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Abstract 

The High-Fidelity Generalized Method of Cells is a new micromechanics model for unidirec- 
tionally reinforced periodic multiphase materials that was developed to overcome the shortcom- 
ings of the original Generalized Method of Cells. The high-fidelity version combines the basic 
elements of the higher-order theory for functionally graded materials with the homogenization 
theory framework, together with a volume discretization approach similar to that employed in 
the original model. The use of a higher-order displacement field approximation in the individ- 
ual subcells of a repeating unit cell representative of a periodic material provides the necessary 
shear coupling absent in the original model, which in turn dramatically improves the accuracy 
of the local stress and strain fields and, in certain cases, the effective moduli. Herein, we test 
the predictive capability of the high-fidelity model in estimating the elastic moduli of periodic 
composites characterized by repeating unit cells obtained by rotation of an infinite square fiber 
array through an angle about the fiber axis. Such repeating unit cells, which typically pos- 
sess no apparent planes of material symmetry, may contain a few or many fibers, depending 
on the rotation angle. In order to analyze such multi-inclusion repeating unit cells efficiently, 
the high-fidelity micromechanics model’s framework is reformulated using the approach previ- 
ously employed by the authors in reformulating the higher-order theory for functionally graded 
materials. The excellent agreement with the corresponding results obtained from the standard 
transformation equations confirms the new model’s predictive capability for a class of periodic 
composites characterized by non-standard, multi-inclusion repeating unit cells lacking planes of 
material symmetry. Comparison of the effective moduli and local stress fields with the corre- 
sponding results obtained from the original Generalized Method of Cells confirms the need for 
the new high-fidelity model, and dramatically highlights the original model’s shortcomings for 
a certain class of unidirectional composites. 
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1 Introduction 


Micromechanical modeling techniques make it possible to efficiently analyze, relative to the actual 
synthesis and experimental characterization of new material designs, the manner in which different 
microstructural details affect the average and local responses of multiphase materials. This is 
important in the development of optimum material architectures for specific applications which 
may require many iterations. Many micromechanical techniques have been developed during the 
past several decades and these include simple Voigt and Reuss hypotheses, self-consistent schemes 
and their generalizations, differential schemes, concentric cylinder models, bounding techniques, 
and approximate or numerical analyses of periodic arrays of inclusions or fibers in the surrounding 
matrix phase (see, for instance, the monographs by Christensen [1], Aboudi [2], Nemat-Nasser and 
Hori [3] , Kalamkarov and Kolpakov [4] , and Markov and Preziosi [5] ) . 

The Generalized Method of Cells (GMC) is an approximate micromechanical model for the 
response of multiphase periodic composites in the elastic and inelastic regimes. It was originally 
developed by Paley and Aboudi [6] using the original Method of Cells proposed by Aboudi [2] as a 
basis, and subsequently extended by Aboudi and others (cf. review article by Aboudi [7]). GMC has 
been made computationally efficient by Pindera and Bednarcyk [8] for unidirectionally reinforced 
composites, and by Bednarcyk and Pindera [9] for composites with three-dimensional reinforcement, 
through reformulations which take advantage of the method’s lack of shear coupling discussed in 
more detail below. These reformulated versions, which provide the enabling technology for efficient 
analysis of periodic composites with detailed (i.e. realistic) microstructural representation, are 
currently employed in NASA’s software called MAC/GMC [10]. 

GMC is based on a geometric model of a rectangular repeating unit cell representative of the 
modeled multiphase material’s periodic microstructure. To mimic the given microstructure, the 
repeating unit cell is discretized into subcells wherein the displacement field is approximated by 
a linear representation based on local coordinates centered at the subcell’s midpoint. An aver- 
aging procedure that involves the enforcement of interfacial displacement and traction continuity 
conditions between the individual subcells applied in a surface-average sense, and periodicity con- 
ditions between adjacent repeating unit cells, establishes the effective or homogenized properties 
of the composite in terms of Hill’s strain concentration factors [11], and geometric and material 
properties of the individual subcells. 

As demonstrated in previous studies on unidirectional metal matrix composites, GMC is typ- 
ically quite accurate at the macro-level under most (but not all) loading situations, but generally 
not very accurate at the micro-level. This is due to the absence of so-called shear coupling which 
provides the required bridge between macroscopically applied normal (shear) stresses and the result- 
ing microscopic shear (normal) stresses necessary for an accurate estimate of micro-level quantities. 
This absence is due to the linear displacement field approximation within each subcell, which pro- 
duces constant strain and stress fields, and also to the manner of applying traction continuity 
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conditions in a surface- average sense at the interfaces separating adjacent subcells. This leads 
to constant traction components along the rows and columns of the repeating unit cell, resulting 
in inaccurate local stress fields. Despite this limitation, GMC and its earlier version have been 
successfully employed to predict the initial yield surfaces and subsequent macroscopic response of 
unidirectional metal matrix composites under different types of loading and with different fiber 
arrays characterized by material planes of symmetry (cf. Pindera and Aboudi [12], Pindera et al. 
[13], and Arnold et al. [14] for earlier investigations, and more recently Iyer et al. [15] and Arnold 
et al. [16]). As demonstrated by Arnold et al. [14] and Iyer et al. [15], the success of GMC in 
predicting the response of metal matrix composites with relatively small constituent elastic moduli 
mismatch is its ability to reproduce with reasonable accuracy the second stress deviator which con- 
trols the evolution of inelastic strains in the matrix phase, despite poor estimates of the individual 
stress components. The computational efficiency of the reformulated GMC makes it an attrac- 
tive tool for stand-alone applications as well as applications which utilize GMC as a subroutine in 
large-scale structural mechanics analyses, provided that the method is employed within its limits 
of applicability and its limitations are recognized as discussed below. 

There are several consequences of the shear-coupling’s absence which must be recognized. First, 
the constant tractions within each row and column of subcells, which result from the application of 
traction continuity at the interfaces separating adjacent subcells in a surface-average sense, dictate 
that the greatest normal or shear traction component in the pertinent direction is directly influenced 
by the most compliant subcell. If the modulus of such subcell happens to be zero, as may occur 
in the presence of porosity, then the normal/shear traction associated with the unit normal along 
the direction of travel will be zero, rendering the row/column containing this subcell ineffective 
in supporting the load. The normal/shear stress re-introduction into such row/column of subcells 
does not occur precisely because of the absence of shear coupling (in contrast to the shear lag-type 
analysis employed to explain the re-introduction of load into a broken fiber away from its ends). 
In particular, the presence of porosity however small in a given plane of subcells ensures that the 
transverse shear stress cannot be supported in the entire plane. Thus for a unidirectional composite 
with a unidirectional porosity in place of a fiber, the transverse shear modulus will be zero according 
to GMC. The same anomaly occurs in the presence of a disbond across an interface separating two 
subcells. The above arguments obviously carry over to subcells with very small moduli relative 
to other subcells. Modeling the response of unidirectional composites with debonded fiber/matrix 
interfaces, discontinuous fiber composites as well as woven composites using GMC requires special 
modifications which mitigate the absence of shear coupling, as recently discussed by Bednarcyk 
and Arnold [17], Pahr and Arnold [18], and Bednarcyk and Arnold [19]. Alternatively, a new but 
more computationally demanding model discussed below, which is also available in MAC/ GMC, 
can be employed to circumvent the above limitations, Bednarcyk and Arnold [10]. 

Further, GMC-based analysis of repeating unit cells that do not have two orthogonal planes of 
material symmetry aligned with the cell’s coordinate system in the investigated plane is not adrnis- 
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sible. This is because GMC cannot accommodate material microstructures leading to anisotropic 
effective properties when the individual phases are orthotropic or isotropic. At most, orthotropic 
effective properties can be predicted which require planes of material symmetry when the individ- 
ual phases are orthotropic/isotropic. This is also a direct result of the absence of shear coupling 
in the formulation. Anisotropic effective properties are predicted only when the individual phases 
are anisotropic. Analysis of repeating unit cells that lack planes of material symmetry erroneously 
produces orthotropic properties and therefore is not admissible. We note that this shortcoming has 
been disregarded by some authors, particularly when applying the method to periodic composites 
with locally irregular microstructures. This is only valid in a very approximate sense for a sufficient 
number of inclusions of similar dimensions within the repeating unit cell which produce isotropic 
effective properties in the plane transverse to the reinforcement direction, Orozco [20]. 

In order to overcome the above shortcomings, a new micromechanics model has been developed 
recently by Aboudi et al. [21,22,23] for the response of multiphase materials with arbitrary periodic 
microstructures, named High-Fidelity Generalized Method of Cells (HFGMC) in part because it 
employs a similar volume discretization as the original GMC. The model’s analytical framework 
is based on the homogenization theory, but the method of solution for the local displacement and 
stress fields within the repeating unit cell utilizes concepts previously employed in constructing the 
higher-order theory for functionally graded materials [24]. The homogenization theory is used to 
construct a higher-order displacement field approximation at the local microstructural level of a 
periodic multiphase material in a consistent fashion, and to derive the governing field equations 
and the boundary conditions that the displacement field must satisfy. The actual solution for the 
local displacement, strain and stress fields follows the above-mentioned higher-order theory. The 
higher-order displacement field approximation at the local level provides the necessary coupling 
between the local normal and shear stress fields and the macroscopically applied loading. This cou- 
pling dramatically improves the accuracy of estimating the local stress fields relative to GMC, as 
previously demonstrated for simple periodic microstructures characterized by orthogonal planes of 
material symmetry through exact analytical and numerical solutions [21,22,23]. It also enables ac- 
curate analysis of unidirectional metal matrix composites in the presence of debonded fiber/matrix 
interfaces in a direct manner, Bednarcyk et al. [25], in contrast with the special modifications in- 
troduced into GMC mentioned earlier [17]. Limited data also has been generated for unidirectional 
composites with locally irregular microstructures using HFGMC, Pindera et al. [26] . It is for such 
composites that the power of HFGMC becomes evident due to the importance of shear coupling 
in the presence of locally irregular microstructures as demonstrated by the above study, further 
highlighting the differences between the two models’ predictive capabilities. 

Given the past experience, it is expected that the shear-coupling effect will play a key role in 
correctly capturing the elastic moduli of unidirectional composites with repeating unit cells that 
lack orthogonal planes of material symmetry due to rotation about the fiber axis, rather than locally 
irregular microstructural arrangement. Therefore, HFGMC is employed herein to determine the 
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effective moduli and local stress fields in unidirectional composites characterized by a repeating unit 
cell with a square fiber array loaded by average stresses that do not coincide with the orthogonal 
planes of material symmetry. Such a situation arises when a unidirectionally reinforced composite 
is rotated by an angle about the fiber axis. The repeating unit cell for a particular rotation angle 
typically may contain a large number of fibers (in contrast with just one for the square array), which 
makes HFGMC computationally expensive due to large number of rectangular subcells required to 
model realistic geometric details (such as circular fibers, for instance). Therefore, in order to enable 
modeling of such multi-inclusion repeating unit cells in sufficient detail, the theoretical framework of 
the HFGMC micromechanical model is reformulated using the local/global stiffness matrix approach 
originally proposed by Bufler [27] for the elastic analysis of isotropic layered media, and extended by 
Pindera [28] to layered anisotropic composites. The same approach has been employed by Bansal 
and Pindera [29] in reformulating the original higher-order theory for functionally graded materials, 
resulting in substantial reduction in the governing system of equations. 

The predictions of the reformulated HFGMC model for the elastic moduli of two types of 
unidirectional composites as a function of the rotation angle are compared with the results obtained 
from the transformation equations for an elastic stiffness matrix of an orthotropic material rotated 
by an angle about the fiber axis. Such rotation produces an elastic stiffness matrix in the rotated 
coordinate system that represents a monoclinic material with just one plane of material symmetry 
(the plane perpendicular to the fiber axis). This provides a critical test on the self-consistency 
of the newly developed HFGMC model previously untested in this manner. Comparison with the 
corresponding GMC predictions are also provided to further justify the development of the new 
model. Towards this end, two types of unidirectional composites with the same inclusion content 
have been employed in this investigation which are characterized by radically different fiber /matrix 
moduli ratios representative of a glass/epoxy system and a porous aluminum matrix. 

2 Theoretical Framework of the High-Fidelity GMC Model 

The theoretical framework of the original version of the High-Fidelity Generalized Method of Cells 
has been described in detail by Aboucli et al. [21,22,23], and thus only a brief synthesis will be 
provided herein in order to make it possible to follow the efficient reformulation’s derivation. The 
high-fidelity model combines concepts from the homogenization theory and the higher-order theory 
for functionally graded materials. The homogenization theory is employed to construct the correct 
form of the displacement field representation in the discretized domains of the repeating unit cell 
which represents the periodic material’s microstructure; to identify the governing field equations 
for the local problem of the repeating unit cell; and to construct appropriate boundary conditions. 
The construction of the displacement field is based on a two-scale expansion of the form: 


rtj(x, y) = «oi(x, y) + <$«h(x, y) + <5 2 it 2i (x, y) + ... 


( 1 ) 
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where x = (xi,X 2 ,x$) are the global or macroscopic coordinates, y = (yi, 'IJ 2 , ]J:i) are the local or 
microscopic coordinates defined with respect to the repeating unit cell, and the different order 
terms characterized by the powers of 6 are y-periodic due to the material’s periodicity. The size of 
the repeating unit cell characterized by the parameter 8 is small relative to the overall material di- 
mensions such that yi = Xi/8. This suggests that a unit displacement at the local scale corresponds 
to a displacement on the global scale of the order of 8. 

The above displacement field representation, together with the relation yi = Xi/8 between the 
two spatial scales, leads to the following strain field decomposition for periodic materials given in 
terms of the average and fluctuating strains £jj(x) and e^(x, y), respectively 


£ij = e ii (x)+4(x, y) + 0{6) (2) 

The average and fluctuating (local) strains are derived from the corresponding displacement com- 
ponents Ui and u’i 


-13 


1 ’ 2 [ d Xj dxi } 


, 1 du- du- 


(3) 


The above strain decomposition makes it possible to express the displacement field in the form 


rti(x, y) = EijXj + u'i + 0(8 2 ) (4) 

where Eij are the known or applied macroscopic strains. This form is employed in constructing an 
approximate displacement field for the solution of the cell problem following the methodology of 
the reformulated higher-order theory for functionally graded materials. 

For specified values of the average strains, the unknown fluctuating displacements must satisfy 
the local stress equilibrium equations subject to periodic boundary conditions imposed on the 
displacement and traction components that are prescribed at the boundaries of the repeating unit 
cell. In addition to these boundary conditions, it is also necessary to impose the continuity of 
displacements and tractions at the internal interfaces between the phases within the repeating unit 
cell. The manner of determining the fluctuating displacements in the repeating unit cell follows the 
methodology employed in the higher-order theory for functionally graded materials. The solution 
is based on satisfying the local stress equilibrium equations in a volumetric sense in the individual 
subvolumes into which the repeating unit cell is discretized in a manner that mimics the actual 
microstructure of the periodic multiphase material. In addition, the interfacial displacement and 
traction continuity conditions are imposed in a surface-average sense at the interfaces separating 
the individual subvolumes, as are the periodic boundary conditions. 

In the original formulation of the high-fidelity theory, a two-level discretization of the repeating 
unit cell’s microstructure was used, involving division into generic cells which were further sub- 
divided into four subcells in the case of periodic materials with continuous reinforcement along 
a common direction. This two-level discretization process, based on the volume discretization 
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employed in the original version of the higher-order theory for functionally graded materials, un- 
necessarily complicated the derivation of the volume-averaged stress equilibrium equations, which 
were obtained in a circuitous manner by satisfying the different moments of the equilibrium equa- 
tions. These manipulations, in turn, suggested that the higher-order theory is a variant of a 
micropolar, micromechanics-based continuum theory, which is not the case. The derivation of the 
traction continuity conditions between adjacent generic cells was also complicated by the two- level 
discretization of the repeating unit cell. 

The reformulation of the higher-order theory was based on a single-step discretization process 
involving only subcells as the fundamental building blocks of a functionally graded material’s mi- 
crostructure. This simplification facilitated the construction of a local stiffness matrix relating the 
surface-averaged displacements to the corresponding surface-averaged tractions of the subcell, with 
the surface-averaged displacements becoming the fundamental unknown quantities. As a result, 
the number of equations governing the micro-level displacement and stress fields within the graded 
material was reduced by as much as 60 % upon assembly of the individual local stiffness matrices 
into the global stiffness matrix [ 29 ] . Furthermore, the reformulation simplified the derivation of the 
governing equations and revealed the higher-order theory to be an approximate elasticity technique 
involving direct volume- averaging of the stress equilibrium equations in conjunction with the impo- 
sition of displacement and traction continuity conditions in a surface- average sense across adjacent 
subcell faces. Herein, the same approach is employed in reformulating the solution for the local 
problem of the repeating unit cell within the High-Fidelity GMC framework. 


3 Efficient Reformulation of the Cell Problem 


The local analysis is performed on the repeating unit cell representative of a periodic material’s 
microstructure in the X2 — x% plane, Fig. 1 (a), and characterized by continuous reinforcement 
along the x\ axis. The periodic microstructure is made up of any number of phases arbitrarily 
distributed within the repeating unit cell so as to produce fully anisotropic behavior in the X2 — x% 
plane. The repeating unit cell is appropriately discretized into subcells, designated by (/?, 7), so 
as to mimic the material’s periodic microstructure, as shown in Fig. 1 (b) for the repeating unit 
cell highlighted in Fig. 1 (a). In this case, 100 x 100 equally dimensioned subcells were used to 
capture the three reinforcement shapes with sufficient fidelity, noting that such refinement would 
have been computationally prohibitive in the original formulation. The indices f3 = l,...,Np and 
7=1, ..., N, y, which span the repeating unit cell dimensions along the IJ 2 and y% axes, identify the 
(/?, 7) subcell in the 1/2 — 2/3 plane. The subcell dimensions along the ij2 and 2/3 axes are hp and Z 7 , 
respectively, such that 

Np iv 7 


n = Y, h e 


and 


L — ^ 


P = 1 


7=1 


where H and L are the overall repeating unit cells dimensions along these axes. 
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Figure 1 . (a) A continuously reinforced multiphase composite with a periodic microstructure 
in the X2 — £3 plane constructed with repeating unit cells, (b) discretization of the repeating unit 
cell into subcells employed in the reformulation of the High-Fidelity GMC. 
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Following the displacement field representation within the repeating unit cell of a periodic 
material given by Eq. (4), the displacement field in each (/?, 7) subcell is written as follows 


1/3, 7) - . / 

U\ = EijXj + II, 


1/3,7) 


(5) 


The fluctuating components of the displacement field arise due to the heterogeneity of the 

medium and are a function of the microscopic (local) coordinates y 1/3,7) = (y^\y^) attached to 
the subcell’s center as shown in Fig. 1(b). These fluctuating components are approximated in 
each (/?, 7) subcell by the same second-order, Legendre-type polynomial expansion in the local 
coordinates as that employed in the original formulation 


' 1/3,7) 
U { 


w(/ 3,7) , y(/3) ixd/3,7) , 4. 

VV U 00) + W2 ^i(lO) + y 3 vv i( 01) + 


*(00) ' ”2 rr i( 10 ) ~ "3 yr i( 01) 

l( 35 f 2 - k l)W^ + 1(35<^ 2 - f )<-> 


i = 1,2,3 


( 6 ) 


where are the unknown microvariables associated with each subcell. Using the above flue- 

lymn) ° 

tuating field representation, the strain components in each (/3, 7) subcell are derived from 


in the form 


p(/3,7) 

ij 


£■ 


h d < 


1/3,7) 


dyf 


+ 


du! 


1/3,7) 


to(-) 


(7) 


,.1/3,7) 

fc ll 


£ll 


£ 

£ 

£ 

£ 

£ 


1/3,7) 

22 

1/3,7) 

33 

1/3,7) 

12 

1/3,7) 

13 

1/3,7) 

23 


£22 + luS + 3yr^ 


2 ( 10 ) 


1/3)ti/(/3,7) 


2(20) 


? _l T/T^l/3,7) 1 QttU) t;i/(/3>7) 
£33 + 1 T 3(01) + 3 y s H/ 3ro2 ^ 


3(02) 


72 + l[< 1 „t + 3# l H' 1 ( & 7 ) 1 ] 
g i3 + 5 Km + 

®23 + + -77 ll 4o/l + '7 < fl0) ) 



( 8 ) 


Using the 


above relations in the Hooke’s law for an orthotropic elastic (/?, 7) subcell 


1/3,7) _ ^1/3, 7) 1/3,7) 
ij ~ ijkl L kl 


(9) 


the subcell stress components are then expressed in terms of the unknown microvariables and 

L 1 i{mn) 

the applied macroscopic strains Eij. For an orthotropic subcell, the stiffness tensor elements 
are characterized by nine independent components in the material principal coordinate system 
formed by the intersections of three orthogonal planes of material symmetry. 

In the original formulation, the unknown microvariables were determined by satisfying 

the local stress equilibrium equations in a volumetric sense, and the displacement and traction 
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continuity conditions between subcells and generic cells, together with the periodic boundary con- 
ditions, in a surface-average sense. In the reformulation, the surface- averaged fluctuating displace- 
ment components become the fundamental unknowns which are related to the surface-averaged 
tractions through a local subcell matrix constructed in the manner described below. 


3.1 Local Stiffness Matrix 


The local stiffness matrix for each (/?, 7 ) subcell is constructed by first expressing the tractions 
defined at the subcell’s outer faces in terms of stresses through Cauchy’s relations 


'i ~ ji n 3 


( 10 ) 


where are components of the unit normal to a given face of the (/?, 7 ) subcell. The surface- 

averaged tractions t^’ 7 ^ are defined at the outer faces of the subcell in the following manner 


f(/3,7) — \p+ P~ P+ P~ P+ P~ l3— l3+ l3— l3+ l3— 1(/3,7) 
l — [i.j , Li ,1*2 1 b 2 ’ l 3 1 L 3 > b l i b \ i l 2 > b 2 ? b 3 > b 3 J 


where 


^± 08 , 7 ) = 1 


y// , n(/3,7) / 1 —(7) \ 1 -(7) 

h (±Y ,y z ) d V3 


^3±08,7) = 


7 J —I7 / 2 

rhp/2 


;= 1,2,3 

% i-V 2 2 

and (j = 2, 3) denotes the direction of the normal to the positive (+) or negative (-) face of the 
(/?, 7 ) subcell, Fig. 2. 


( 11 ) 

( 12 ) 

(13) 



Figure 2. A view of a subcell illustrating the convention employed in designating the surface- 
averaged displacement and traction components employed in the reformulation. 
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Substituting Eqs. (8) — (10) into Eqs. (12) — (13) , performing the required integration and 
assembling the resulting equations in matrix form, we obtain 

1 (&7) 


i3+ 
t 3 
f 3- 
z 3 


1 (0.7) 


_ (0.7) 

°33 


+2+ 

Z 1 

(0, 7) 

E; 

7a 

11 

1 

3hp 

2 


W 1(10 ) 

f 2 - 

L '1 J 


-1 

3 hp 
2 


^1(20) 

i3+ 

'] 

(0,7) 

TCI lO 

X? 

11 

1 

3l~/ 

2 


VEi(oi) 

+3“ 

z i 


-1 

3 Z 7 
2 


^1(02) 

X2+ 

l 2 

(0,7) 

7a cm 

0 " 

11 

1 

3h fJ 

2 


E^io) 

t 2 - 

z 2 


-1 

3hp 

2 


^2(20) 

J2+ 

z 3 

(0,7) 

7 a ^ 

11 

1 

3/1/3 

2 


^3(10) 

f 2 - 

L Z 3 J 


-1 

3hp 

2 


^3(20) 

i3+ 

l 2 

(0,7) 

7 a ^ 

11 

1 

3i 7 

2 


^2(01) 

z 2 


-1 

3l~f 

2 . 


^2(02) 


I o/^<(0>7) 

“r -^66 


+ 2C<f 1 1 


£12 
— £12 

£13 

£13 


(0.7) 


■G 


(0,7) 


23 


IE 


3(01) 


-IE 


3(01) 


(0.7) 




(0,7) 

i2 


i= 1 


1 (0,7) 


I /^(0,7) 
~r O44 


TE 


1 (0,7) 


2 ( 01 ) 


—IE 


+ 2Cif 7) 


i /-«(0,7) 
+ u 44 


IE 


2(01) J 

1 (0,7) 


3(10) 

|_ —^3(10) J 


+ 2C|f 7) 


(14) 

(15) 

£ii 

(16) 

£23 
-£23 

(17) 

£23 
“£23 

(18) 


, 3/ 7 

1 2 


^3(01) 

1 3i 7 

1 2 


^3(02) 


1 (0,7) 


+ c<E 


w, 


1 (0,7) 


2 ( 10 ) 


-IE 


2(10) J 




(0,7) 

i3 


i=\ 


Si 

-Si 


(19) 


Equations (14) — (19) relate the surface- aver aged tractions to the first- and second-order mi- 
crovariables. Next, we express the microvariables E 7 ^ 7 ^ in terms of the fluctuating surface- 
averaged displacements u'^ 3,7 ) 


u 


,(3 ■v) r — / 2— |— _ / 2 — _ / 2+ _ / 2 — — / 2— |— _/2 — — / 3— |— _/3 — _ / 3-|- _/3 — _ / 3-b _/3 — if/5'7) 

— 71 'll n, n, ?/„ 7/ 7 / 7 /_ ?/_ n, 01 VA-'W/ 


( 20 ) 


where 




u. 


2±(0, 7 ) 


3±(0,7) 


= I [ h/ * u'V"\±^A' l) )dvi' l) 


J—l-y/2 


rhg/2 


h. 


u. 


(±f .i/rw 


(21) 

( 22 ) 


i = 1,2,3 

10 4-^/2 

and where u( J± ^’ 7 ^ is the fluctuating surface-averaged displacement in the i th direction evaluated on 
the face of the (/?, 7) subcell with normal in the ±j th direction, Fig. 2. Substituting the expressions 
for given by Eqs. (6) in the above definitions, performing the averaging procedure, and 

assembling the resulting fluctuating surface-averaged displacements in matrix form yields 


r - /2 + 1 

( 0 , 7 ) 

«i 


-, 2 - 

V u i \ 



hp h% 


2 4 

hp 

2 4 


lEi 


1 (0,7) 


1 ( 10 ) 


IE 


1(20) J 


W^ 


1 ( 0 , 7 ) 


1 ( 00 ) 


[ ^1(00) J 


( 23 ) 
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— /3+ 

(Pa) 

u i 


_/3— 

L u i \ 


r -' 2 + i 

(Pa) 

U 2 


Si 

10 to 


r -' 2 + i 

u 3 

(Pa) 

- 1 2 - 

l U 3 J 


— /3+ 

(Pa) 

u 2 


_/3— 
_ u 2 


-/3+ 

u 3 

(Pa) 

-,3— 

L U 3 J 



hp_ 

2 

hf 3 

2 

hp_ 

2 

hp_ 

2 


r 2 
lx 
4 

i 2 

IX 

4 


/i 2 

n 0 

4 

h 2 

n p 


h 2 

n i3 

4 

h 2 

n l3 

4 


l 2 

lx 

4 

l 2 

lx 

4 


; 2 

lx 

4 

Z 2 

lx 

4 


Wi 


(Pa) 


1 ( 01 ) 


Wi 


+ 


Wi 


1 ( 00 ) 


1(02) 


W, 


1 ( 00 ) 



1^2(10) 

(Pa) 

+ 

w 2m 


W / 2(20) 

H2(00) 


ll3(io) 

(Pa) 

+ 

113(00) 


H / 3(20) 

113(00) 



H2(01) 

(Pa) 

+ 

1^2(00) 


ll / 2(02) 

H^2(00) 


113(01) 

(Pa) 

+ 

113(00) 


113(02) 

H3(00) 


(P, 7) 


(0,7) 


1 (0,7) 


(Ax) 


(Ax) 


Adding and subtracting rows in each of the above equations, we obtain the following 

i (Pa) 


to 

0 0 

(Pa) 

cn' 

0 0 

£ £ 

(Pa) 

112(10) 

112(20) 

(Pa) 

113(10) 

H3(20) 

(Pa) 

112(01) 

112(02) 

(Pa) 

H3(01) 

113(02) 

(Pa) 


hf 3 
2 

ft? 


l 7 

2 

zj 


L 

j_ 

hp 

2 

ft? 


l 7 

2 

Zj 


2 

7T 

l 7 


Zl/3 

2 


2 

Zj 



r -' 2 + 1 

u i 

(Pa) 

4 

0 


-1 2 - 

L u i J 

n p 

Hl(OO) 



— /3+ 
«i 

(Pa) 

4 

0 


_/3— 

L u i J 

~ Z 2 

Hi (oo) 


(Ax) 


’ hp 
2 


2 

Zj 


2 

F 



r - /2 + 1 

u 2 

(Pa) 

4 

0 


SSl 

M i2 

h p 

11 2 ( 00 ) 


.-4 + 1 

(Pa) 

4 

0 


Si 

1 

_ h2 
n 0 

H 3 (oo) 


(Ax) 


1 (Ax) 


— /3+ 

n 2 

_/3— 

Ur, 


_/3+ 

u 3 

_/3— 

Un 


(Pa) 


(Pa) 


1 2 
4 X 


/2 

s 




2(00) 


(Pa) 


(Pa) 


W, 


(24) 

(25) 

(26) 

(27) 

(28) 
relations 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


3(00) 


Equations (29) — (34) relate the first- and second-order microvariables to the fluctuating surface- 
averaged displacements and the zeroth-order microvariables. In order to express the microvariables 
explicitly in terms of the surface-averaged quantities, the local stress equilibrium equations are 
satisfied in a volumetric sense. We note that the reformulation eliminates the need to consider the 
first and second moments of these equations, simplifying the volume averaging procedure of the 
equilibrium equations to the zeroth moments as shown below 
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Expressing the subcell stresses in terms of the fluctuating surface-averaged displacement using 
the expressions for the first- and second-order microvariables given by Eqs. (29) — (34) in 

Eqs. (8) and the resulting relations in Eqs. (9), and performing the required integration in the 
above equation yields after simplification 


where 


= ^kr(“. 2+ + 7~) ,M + + «i~) w 


(0,7) 


h 2(i(0,7) 

llaLscc: 


1 ( 00 ) 


2 C{f’ 7) 


Wr 


(0,7) 


c. 


(0,7) 


2 ( 00 ) 


22_ 

(0,7) ' 
''22 

i(0,7) 


2 Z2ci?’ 7) 

u2 (-1(0,7) 


7) 


(u? + + U2 2_ ) (/3,7) + 13 7 «2 3+ + %”) (/3 ’ 7) 

0 ^, 1 ) v 2 2 ' onn(0,7) v 2 2 7 


— *7-3 /-'3+ -'3-x(/3, 7 ) t 7 L/ 44 /- -/2+ ,2 -x(/ 3 j7 ) 

”3(00) n/=i(0,7) '3 + U 3 1 + 9h 2r(0’7) [ 3 + 3 j 

ZUqq Zr / L qq 


33 


]2ri(0,7) 

>2f,(0,' 

qW 33 


/7(d,7) _ fi(0,7) , ^0_fi(0,7) 

°11 — L '66 + ^2 C '55 


c. 


(0,7) 

22 


(-,( 0 , 7 ) 1 

— <-'22 ■+■ |2 C 44 


7 

hi 


0n(0,7) 
2 
7 


(7(0,7) _ ri(0,7) , 7 /~f(d,7) 

~ L 33 + i,2 <-'44 

n 0 


(36) 

(37) 

(38) 

(39) 

(40) 

(41) 


Substituting the zeroth-order microvariables back into Eqs. (29) — (34) and then the resulting 

equations for the first- and second-order microvariables !e 7’ 7 \ into the equations for surface- 
averaged tractions, Eqs. (14) — (19), we obtain the expressions for the surface-averaged tractions 
exclusively in terms of the fluctuating surface-averaged displacements and the applied macroscopic 
strains. We observe that the inplane surface averaged tractions (i,j = 2, 3) are solely functions 
of the inplane fluctuating surface-averaged displacements ( i,j = 2, 3), and are related through 

the local stiffness matrix as shown below 


12+ 

l 2 

(0,7) 

' K u 

Kn 

0 

0 

Kn 

K u 

77,7 

Iu 8 ' 

(0,7) 

r ~/2+ 
«2 

f 2 - 

l 2 


K 21 

K22 

0 

0 

77 25 

K26 

K27 

7728 


_/2— 

«2 

J2+ 

l 3 


0 

0 

K33 

K34 

K35 

K 3 q 

K 3 7 

7738 


_/2+ 

“3 

l 3 


0 

0 

K43 

K 44 

K45 

K 4 G 

K47 

7748 


_/2— 

«3 

J3+ 

l 2 


k 51 

A'52 

K 53 

7754 

K55 

77,56 

0 

0 


_/3+ 

«2 

l 2 


Kqi 

Kq 2 

Kq3 

Kq4 

Kq5 

7766 

0 

0 


_/3— 

« 2 

J3+ 

73 


K n 

K72 

K73 

K74 

0 

0 

K77 

7778 


_/3+ 

% 

1 

ftS 5 


. K 8 i 

K 8 2 

K 8 3 

K 8 4 

0 

0 

K 8 7 

7788 . 


_/3— 
L «3 
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-I (0,7) 


C\2 C 22 C 23 

— C 12 —C 22 —C 23 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

Cl. 3 C 23 C 33 

— Ci3 — C 23 — C 33 


0 

0 

20*4 

-2044 

2044 

-2044 

0 

0 


£11 

£22 

£03 

£23 


( 42 ) 


where the elements of the local stiffness matrix Kij are given explicitly in terms of the mechanical 
properties and subcell dimensions in the Appendix. Similarly, the out-of-plane surface- averaged 
tractions (j = 2, 3) are related to the out-of-plane displacements u A 3± ( j = 2,3) through the 
local stiffness matrix as shown below 


" *T ' 

(0,7) 

An 

Ll2 

A13 

A14 

(0,7) 

— / 2+ 
u i 

(0, 7) 

Cm 

0 



A21 

^22 

L23 

A24 


_,2— 
u 1 

+ 2 

Cqq 

0 

i3+ 

Z 1 


A31 

L32 

L33 

A34 


_,3+ 

U 1 

0 

C55 

f 3 - 

L J 


L41 

A42 

A43 

A44 


-i 3 — 

l U 1 J 


0 

— C55 


1 (0,7) 


£12 

£13 


(43) 


where the elements of the local stiffness matrix Lij are also given explicitly in terms of the mechan- 
ical properties and subcell dimensions in the Appendix. 


3.2 Global Stiffness Matrix 

The local stiffness matrices are used to construct the global stiffness matrix using first the interfa- 
cial displacement and tractions continuity conditions, and then the periodic boundary conditions, 
all imposed in a surface-averaged sense across adjacent subcell interfaces. We start with the dis- 
placement continuity conditions. Considering the (3 th interface, the fluctuating surface- averaged 
displacements in the x\ direction, u'i + ^ a \ u ) 2 ' ;3+1,7 \ and the corresponding displacements in 
the X 2 and X 3 directions, u '2 and u 3 + ^’ J \ u 3 ^ +1,7 \ respectively, must be 

equal. Hence, we represent them using common unknown quantities, i.e., 


,2+09,7) 

U 1 

,2—09+1,7) 
— «1 

,209+1,7) 

— 

(44) 

,2 +(0, 7 ) 
u 2 

,2—09+1,7) 
— « 2 

,2(0+l,7) 
— « 2 

(45) 

,2+09,7) 

“3 

,2—09+1,7) 

— 0.3 

,209+1,7) 
— u 3 

(46) 
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Similarly, applying displacement continuity at the 7 th interface in the x\, X 2 and X 3 directions, 
respectively, in a surface- average sense we have 


;'3 + (0,7) 

1 

,3-09,7+1) 
— u 1 

,309,7+1) 
— u. 1 

(47) 

/3+(/3,7) 

'2 

-;'3— (d,7+l) 
« 2 

,309,7+1) 
— «2 

(48) 

i'3+(/3,7) 

'3 

-'3— (/3,7+l) 
— “3 

,309,7+1) 
— u 3 

(49) 


Equations (44) — (46) and (47) — (49) hold true for (3 = 1, ...,Np—l and 7 = 1 , _ZV 7 — 1 , respectively. 
This gives rise to 3 (Np — 1 )_ZV 7 + 3 (1V 7 — 1 )Np unknown surface- averaged displacements defined at 
the subcell interfaces along the X 2 and 23 directions. The remaining 6(Np + 1V 7 ) surface-averaged 
displacements at the external boundaries of the repeating unit cell 


,2(1 jT ) 

,3(^,1) 

, 


,2(1, 7 ) ,2(1,7) -/2(AT / 3 +l,7) ,2(^+1, 7 ) ,2(^+1, 7 ) 

a 2 > u 3 ) “1 j a 2 ? u 3 1 

,303,1) -,309,1) ,3(/3,JV 7 +l) ,3(/3,lV 7 +l) ,3(/3,TV 7 +l) 

« 2 , “3 j u l 1 u 2 ’ a 3 


(50) 


are related to the corresponding surface-averaged boundary tractions. These surface-averaged 
boundary quantities are related through the periodic boundary conditions which will be incor- 
porated into the global stiffness matrix. 

Next, the traction continuity conditions at the interfaces between adjacent subcells are imposed. 
Considering the (3 th interface, the continuity of tractions in the x\, X 2 and X 3 directions, respectively, 
is ensured by 


+ f~(P+hl) = Q 

f+(0,7) + f-(0 +hl ) = Q 

+ j 2 -(/Hl, 7 ) = o 


(51) 

(52) 

(53) 


where (i = 1,2,3) are the surface-averaged traction components on the right face (T 2 = 

h.p/2) of the (/?, 7 ) subcell and if ^ +1,7 - > (i = i ) 2, 3) are the surface-averaged traction components 
on the left face (X 2 = — hp + \/2 ) of the {(3 + 1 , 7 ) subcell , defined by Eq. (12). Similarly, applying 
traction continuity at the q t/! interface in the aq, X 2 and X 3 directions, respectively, in an average 
sense we obtain 


£-3+(/3,7) _|_ ^3— (/3,7+l) _ q 

^3+(/3,7) _j_ (9,7+1) _ q 

g+K/3,7) _j_ ^3-(/3,7+l) _ q 


(54) 

(55) 

(56) 


Equations (52)-(53) and Eqs. (55)-(56) can be written in terms of the common interfacial surface- 
averaged displacements u'^' ’ \ ’ \ u : jf' ’ \ u 3 ^ ’ ^ using the local stiffness matrix that relates the 

inplane surface-averaged tractions and displacements, Eq. (42), and the displacement continuity 
conditions given by Eqs. (45)-(46), and Eqs. (48)-(49) 
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K g, 7 )_, 2( / 3 , 7 ) 


(9,7) 


AAA 


+ (^+1 
3(9,7+!) 


+ A' 


09+1, 7 )x -,209+1,7) 


22 


m. 


+ A 


(9+1,7)-, 209+2, 7) 


21 


U. 


+ K. 


26 


+ A 


25 


+ A' 
f A 


(A7),-'3(/3,7) 


16 


18 


U, 


Un 


+ 

'+ 


AT 


(/3.7)^' 3 (/3>7+ 1 ) 


17 


U 


+ K, 


(9+1,7) ,3(9+1,’ 7) 


28 


u 


+ A. 


(9+1,7) ,3(9+l,7+l) _ 

27 “3 


(AA’ - cg +1 ' 7 >)5n + (c£ 7) - cg +1 ' 7l )522 + (AA* - C<^>ye 33 (57) 


AAA 20 ’' 7 ’ + (if. 


(9,7) 

33 


I to(9+1,7)' 

-I" ^44 


_,2(9+l, 7) 


It. 


+ K, 


(9+1, 7)-, 2(9+2, 7) 


43 


U- 


+ A 


(/3,7) J7 '3(9,7) 


36 


U. 


+ 


A 


(9+1,7) ,3(9+1, 7) 

46 “2 


+ a: 


(9,7) 47' 3(9, 7+1) 


35 




+ A, 


(9+1,7)-, 3(9+1, 7+1) 


45 


It. 


+ K^u' mi) 


to( 9+1,7) ,-7'3(9+l>' 7) , 7^(9, 7)„- 7 '3(9, 7+1) ^(9+1, 7) ,3(9+1, 7+1) 

-“-48 ttq “r -“.^7 u -3 ~r 7V47 (+ 


l 37 


L 47 


= 2(C, 


(9,7) 

44 


38 “3 + 

(9+1,7) 


a 


44 


)e 23 (58) 


1^(9, 7) 47' 2(9,7) , ^(9,7)42/2(9+1,7) (9,7+1) ,2(9, 7+1) (9,7+1) ,2(9+1, 7+1) , 

^52 u-9 I -iVj^i 0,9 1 TV ^9 L69 I Tv^-| ^9 I 


51 


L 62 


61 


A, 


(/3,7)„-,'2(9,7) 


54 


U. 


+ A, 


(9,7)47,2(9+1,7) 


53 


u 


+ A, 


(9,7+1) 47, 2(9, 7+1) 


64 


U 


+ A, 


(9,7+1)-, 2(9+1, 7+1) 


63 


It 


+ 


<A7) + (A?’ 7 ’ + A6 ,7+1) )“2 3|,5 ’ 7+11 + Af 7+1) »2 3<A7+2) = 2(C<J' 7) - C<f 7+1) )£23 (59) 


AAA 20J ' 7) + AAA 2<,S+1 ' 7) + AA +1 A 2(ft7+1) + AA+'A^-^+A 


A 


L 72 

(9, 7)47/2(9,7) 


74 


U 


+ a: 


L 71 

(9,7)47/2(9+1,7) 


81 


73 


U 


^ 78 , 7 ) “ 3 3(/3 ’ 7) + ( A " 


77 


+ -^84 

+ a; 


(9,7+1)-, 2(9, 7+1) 


17 


(9, 7+1), 77' 2(9+1, 7+1) 


Ito 


+ ^83 

+ a 87 


17 


+ 


Uo 


,,o(9,7) ^(9,7+1) \= 1 t/o(/?,7) -o(9, 7+1) \ ? , />o(9, 7) /o(9,7+l)\ = 

A 13 A 13 All + (0 23 G 23 j£22 + (b 33 C- 33 A 33 


(60) 


Similarly, Equations ( 51 ) and ( 54 ) can be written in terms of the common interfacial surface- 
averaged displacements it) , rt^ ’ J using the local stiffness matrix that relates the out-of-plane 
surface-averaged tractions and displacements, Eq. ( 43 ), and the displacement continuity conditions 
given by Eqs. ( 44 ) and ( 47 ), 


L (9,7) -/2 ( 9,7) + (L (9,7) 


J 11 


+ A 


(9+1,7) ^ _ , 2(9+ 1 ,7) ^(9+1,7)47,2(9+2,7) 


22 


21 


U , 


+ A 


(9,7)47/3(9,7) 


14 


It 


+ 


7 (9’7) 47/ 3(9,7+!) , r (9+l,7)47'3(9+l,7) , t (9+1,7) 47/3(9+1,' 7+1) _ 0,70(9,7) 70(9+1, 7)f ? ,4-1 \ 

A13 +A 24 U 1 +A 23 % - A c ' 66 -A 66 )£12 (61) 


r ( 9 , 7)47/2(9,7) 

“32 U 1 


+ A 


(9,7) _, 2(9+1, 7) ^(9,7+1)47,2(9,7+1) _j_ r ( 9 , 7+1)47/2(9+1, 7+1) 


31 


U 


42 


U 


r( 9 , 7 ) 47 / 3 ( 9 , 7 ) 
“34 U 1 


(A. 


(9,7) 

33 


I 7- (9,7+1) 5 
+ -“44 ) 


,3(9, 7+1) + ^(9, 7+1)47/3(9, 7+2) 


43 


It, 


+ L 

= 2(G 


41 

(9,7) 

55 


' + 


c, 


(9,7+1)57= 


55 


)ei 3 ( 62 ) 


Equations ( 57 ) — ( 60 ) and ( 61 ) — ( 62 ) provide us with a total of 3 (Np — 1 )_/V 7 + 3 ( 1 V 7 — l)JV)j equa- 
tions in terms of the common interfacial surface-averaged displacements and the surface-averaged 
displacements at the external boundaries of the repeating unit cell. 


Assembly and structure of the global stiffness matrix We first assemble the global stiffness 
matrices by assuming that the repeating unit cell’s boundary is subjected to prescribed surface- 
averaged tractions. The assembly thus includes the boundary cells (1,7), (Np, 7) and (/?, 1 ), (/?, 1 V 7 ), 
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which provide the additional 6 (Np + N- f ) equations involving the boundary surface- averaged trac- 
tions and displacements, in addition to the application of the interfacial continuity conditions 
described above for the inplane and out-of-plane surface-averaged interfacial displacements, Eqs. 
(57)-(60) and Eqs. (61)-(62), respectively. The final system of equations relating the inplane 
quantities is symbolically written as follows 

' Kn 0 K 13 K 14 1 r ul, 2 1 r t| 1 r ACn AC 12 AC 13 0 1 r en ' 

0 K 22 K 23 K 2 4 Ug 2 0 0 0 AC 24 £22 

. q — _ 0 ~\r (Oo ) 

K 3 i K 32 K 33 0 u 2 t 2 0 0 0 AC 34 e 33 

K 44 K 4 2 0 K 44 Ug 3 tg AC 4 i AC 42 AC 43 0 £23 

where 


_/ 2 

u 2 = 

[u' 2(11 ,. 

_/ 2(2V-y)n 
,.,U 2 J 

with 

of 1 ’ = 

K 2,1 ’ 7) ,.. 

/2(A^3+1,7). 
.., U 2 

_/2 

Ug = 

[uf 11 .. 

;2(W 7 )i 
• • > u 3 J 

with 

of 7 ’ = 

[of’ 71 ,.. 

/2(AT /3 + 1 ,7). 
.., U 3 

_/3 

u 2 = 

[u' 3(11 ,. 

,3(iV /3 ) 1 

■•) u 2 J 

with 

= 


/3(/3,AT 7 +l) 

_/3 

Ug = 

[uf 11 ,. 

,-,'3(^1, 
u 3 J 

with 

u?'» = 


,309,JV 7 +1) 

.., U 3 


and the structure of the surface- averaged traction vectors t 2 , tg, t 2 , t| is similar to the above 
surface-averaged interfacial displacement vectors. In this case, however, the only non-zero surface- 
averaged traction components are those associated with the external surfaces of the boundary 
subcells as shown below. This follows from the interfacial traction continuity conditions given by 
Eqs. (52)-(53) and Eqs. (55)-(56) in the case of the inplane quantities. 


F 2 - 
l 2 ~ 


, 2(W 7 )i 

"A 2 J 

with 

* 2 ( 7 ) 

x 2 

II 

0 

T-2(iV0+l,7)- 
,.,u, t 2 

F 2 - 

t 3 — 


r2(AT 7 ) 1 

■•A3 J 

with 

r2( 7 ) 

z 3 

= i*1 (1 - 7, ,o... 

T-2(iVa+l,7)- 

■) u A 3 

F 3 - 

t 2 — 


tHNp), 

"2 J 

with 

r3(/3) 

l 2 


n f3(/3,W 7 -|-l) 
. ., U, f- 2 

rM 

coco 

II 


^(ATa) 

■•A3 J 

with 

*308) 

z 3 

II 

oTclo 

5 

0 

n f3(/3,W 7 +l) 

. u , (g 


The global stiffness matrix in the above system of equations consists of twelve non-zero subma- 
trices. The diagonal submatrices Kn, K 22 , K 33 and K 44 relate the surface-averaged tractions to 
displacements along their respective directions (for instance Kn relates t\ to u 2 , K 22 relates t| to 
Ug 2 , etc..) and have entries concentrated along the diagonal. The off-diagonal submatrices couple 
the surface-averaged quantities along the x 2 and x 3 directions (K 43 relates t| to u 2 3 , K 44 relates t 2 
to Ug , etc..) and have entries scattered throughout. The structure of the diagonal and off-diagonal 
submatrices, viz. Kn, K 43 etc. constituting the inplane global stiffness matrix has been described 
in detail by Bansal and Pindera [29] in the context of reformulated higher-order theory. The size 
of the global stiffness matrix for the inplane surface-averaged interfacial displacements is 

[21V1V 7 + 1) + 2 N y (Np + 1)] x [2Np(N 1 + 1) + 2N 1 (N p + 1)] 
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Finally, the column submatrices ACn, , , , AC44, which are multiplied by the macroscopic strains 
on the right hand side of Eqs. (63), represent the differences in the elastic stiffness elements 
between adjacent subcells in the 1/2 and 2/3 directions, as shown on the right hand side of the traction 
continuity equations given by Eqs. (57)-(60). 

Similarly, the final system of equations relating the out-of-plane quantities is symbolically writ- 
ten as 
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and the structure of the surface-averaged tractions t 2 and t 3 is similar to the corresponding inplane 
quantities. The global stiffness matrix in the above system of equations consists of four submatrices. 
The diagonal submatrices Ln, L22 relate the surface-averaged tractions to displacements in their 
respective directions i.e. , t 2 to u ^ 2 and t 3 to u ) 3 , respectively, and have entries concentrated along 
the diagonal. The submatrices L12, L21 represent coupling of the surface-averaged quantities in 
the X2 and X3 directions, t 2 with u ) 3 and t 3 with u) 2 , respectively, and have entries scattered 
throughout. The structure of the submatrices Ln and L22 is similar to the structure of Kn and 
the structure of L12 and L21 is similar to the structure of K43. The size of the global stiffness 
matrix for the out-of-plane surface-averaged interfacial displacements is 


[^(^7 + !) + ^7(^0 + 1)] X + 1) + ^7(^/3 + !)] 

Finally, the column submatrices Acn, , , , AC22, which are multiplied by the macroscopic strains 
on the right hand side of Eqs. (64), represent the differences in the elastic stiffness elements 
between adjacent subcells in the 3/2 and 1/3 directions, as shown on the right hand side of the traction 
continuity equations given by Eqs. (61)-(62). 


Reduction of the global stiffness matrices The global stiffness matrices given in Eqs. (63) 
and (64) are further reduced using periodicity conditions on the surface-averaged displacements 
and tractions imposed on the external surfaces of the boundary subcells around the repeating unit 
cell. The periodicity conditions for the surface-averaged boundary displacements are 
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Similarly, the periodicity conditions for the surface-averaged boundary tractions are 
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(65) 


( 66 ) 
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The imposed periodicity conditions (65) — (66) eliminate the traction vectors on the right hand 
sides of Eqs. (63) and (64), and provide us with the necessary 3 (Np + 1)_/V 7 + 3(IV 7 + 1 )Np 
relations for the 3(Np + l)iV 7 + 3(_/V 7 + l)Np unknown subcell surface- averaged displacements, i.e. , 
3(Np — 1)_/V 7 +3(iV 7 — l)Np unknown common interfacial surface-averaged displacements along with 
G(Np+N~/) unknown surface-averaged displacements at the external boundaries of the repeating unit 
cell. These relations are obtained from Eqs. (63) and (64) by combining and deleting appropriate 
rows and columns of the original stiffness matrices appearing in these equations. 

The resulting reduced equations relate the unknown surface-averaged interfacial and boundary 
displacements to the applied macroscopic strains through the reduced stiffness matrix. The final 
reduction of this singular matrix involves constraining the corner subcell faces to eliminate rigid 
body motion. In view of the imposed periodicity conditions on the surface-averaged boundary 
displacements, constraining the external surfaces of one corner subcell and just one appropriate 
external surface of two corner subcells at opposite ends of the diagonal is sufficient. 


3.3 Homogenized Constitutive Equations 


The average strains in each subcell are related to the average macroscopic strains through the strain 
concentration tensor, Aboudi et al. [21], 


^(/3,7) = g 


(67) 


The average subcell strains are obtained by averaging Eqs. (8) over the subcell volume, yielding 
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( 68 ) 


The solution of the reduced systems of equations yields the interfacial and boundary surface- 
averaged displacements as a function of the macroscopic strains. We then calculate the microvari- 
ables associated with each subcell in terms of the macroscopic strains using Eqs. (36) — (38) and 
Eqs. (29) — (34) . This makes possible to obtain the average strain components in each subcell in 
terms of the macroscopic strains from Eqs. (68). 

In order to determine the elements of the strain concentration tensor A^’ 7 ) for each subcell, 
only one component of the macroscopic strain e is applied at a time when solving the reduced 
systems of equations. For instance, applying £ n = 1 and all others zero, and then solving the 
reduced systems of equations to obtain e^d,7) f or each subcell, we obtain the first column of the 
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strain concentration tensor using Eq. (67) . By successively applying the remaining macroscopic 
strain components one non-zero component at a time, we obtain the remaining elements of the 
strain concentration tensor for each subcell. 

The average stress in each subcell is given by 

„(f3n) = c (/3,7)^(/3,7) (69) 

in view of the fact that the material occupying a given (/5, 7 ) subcell is homogeneous. Using Eq. 
(67) to express in terms of the macroscopic strains, we obtain 

7) = (70) 


Averaging the subcell stresses over the entire repeating unit cell, we then obtain the macroscopic 
stress in terms of the macroscopic strains for the composite in the form 
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(71) 


which can be written as 


a = C* e 


(72) 


where C* represents the effective elastic stiffness matrix for the repeating unit cell and is given by 


C* 
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(73) 


4 Numerical Results 

We test the high-fidelity model’s predictive capability by determining the effective moduli of a 
unidirectional composite, with a square array of fibers in the X 2 — X 3 plane, as a function of 
the rotation angle 6 about the fiber axis xq. The moduli in the rotated coordinate system are 
then compared with the standard transformation equations which provide the correct answer. We 
also generate the local stress fields within the repeating unit cell for the given rotation angle. To 
highlight the advantages and need for HFGMC, the moduli as well as the local stress fields predicted 
by this model are compared with the corresponding GMC results. This comparison illustrates the 
importance of including the effects of shear coupling in heterogeneous materials. 

Figure 3 shows the investigated square array of fibers with a fiber volume fraction of 0.35, 
extending to infinity in the X 2 — X 3 plane. Both the fiber and the matrix phases are isotropic. 
To amplify the influence of shear coupling, we consider two cases with radically different contrast 
between the fiber and matrix properties. In the first case, the matrix is an epoxy resin and the 
fibers are glass with typical elastic moduli that produce the Young’s moduli ratio Ef/E m = 20. In 
the second case, we consider an aluminum matrix weakened by holes which are simulated by very 
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Figure 3. A representation of an infinite, square-packed array of inclusions, showing five repeating 
unit cells of the same array in different coordinate systems rotated by the indicated angles about 
the fiber axis, arranged according to increasing number of inclusions contained within each cell. 


Material 

Young’s Modulus (MPa) 

Poisson’s ratio 

Glass fiber 

70, 000 

0.25 

Epoxy matrix 

3,500 

0.35 

Compliant fiber 

700 

0.33 

Aluminum matrix 

70, 000 

0.33 


Table 1. Material properties of the fiber and matrix constituents. 
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compliant inclusions that yield the Young’s moduli ratio Ef/E m = 0.01. The actual constituent 
moduli values are given in Table 1. We note that the Young’s moduli ratios for ceramic and metal 
matrix composites fall well within these limits, and therefore the effective moduli of these materials 
will not be as severely affected by lack of shear coupling. 

As shown in Fig. 3, five different repeating unit cells are analyzed which produce homogenized 
properties of the same fiber array relative to five coordinate systems generated by rotating the 
principal material coordinate system through an angle 8 about the fiber axis. These are arranged 
in two rows such that the number of fibers in each repeating unit cell increases in each row from left 
to right. As observed, the rotation angle does not increase monotonically with increasing number of 
fibers. The first repeating unit cell in the first row of Fig. 3 with the circular fiber in the center, for 
which 8 = 0 °, represents the infinite square fiber array in the principal material coordinate system 
in light of the fact that the orthogonal planes of material symmetry passing through the fiber center 
coincide with the global x\ — X 2 — £3 coordinate system. Clearly, this is the simplest repeating 
unit cell. In order to simulate the response of this square fiber array in the coordinate system 
rotated by 14.04° about the x\ axis, the first repeating unit cell shown in the second row of Fig. 3 
containing a total of 17 fibers is employed. The repeating unit cell representing the composite in 
the coordinate system rotated by 26.57° with respect to the principal material coordinate system 
is the third cell in the first row of Fig. 3. In this case, there are just 5 fibers. The repeating unit 
cell for the rotation angle of 36.87° that appears last in the second row of Fig. 3 contains 25 fibers, 
while the one for the rotation angle of 45.0° which appears second in the first row of Fig. 3 contains 
just 2 fibers. These repeating unit cells were constructed by connecting the center of a reference 
fiber to the center of the fiber a certain number of fibers to the right of the reference fiber and up. 
This specified both the rotation angle and the length of the lower inclined edge of the repeating 
unit cell. Completing the square in the same manner produced the entire repeating unit cell for 
the particular rotation (or inclination) angle. Thus the four rotation angles were obtained from the 
relations 


8 


_i r l 1 3 1, 
tan 1 — , — , — , — 1 
V 2’ 4’ I s 


where the denominator represents the number of fiber distances to the right of the reference fiber, 
and the numerator the number of fibers up. The fiber distance is the horizontal or vertical distance 
between adjacent fiber centers. It can easily be observed that the generated repeating unit cells are 
the basic building blocks of the same fiber array in the five considered coordinate systems, which 
include the principal material system. It is also clear that three out of the five do not possess planes 
of material symmetry. This will result in anisotropic behavior in the X 2 — X 3 plane, necessitating the 
use of periodic boundary conditions easily accommodated by the high-fidelity model’s framework. 

The actual repeating unit cells in the four rotated coordinate systems used in the calculations are 
shown in Fig. 4 in discretized form. Table 2 provides information on the actual microstructural dis- 
cretization used for each repeating unit cell, and the volume fraction of the fiber phase. We note that 
the fiber fraction of each repeating unit cell varied slightly from the nominal fraction of 0.35 due to 
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(a) 6 = 45° (1 00 x 1 00 subcells) 



0 0.2 0.4 0.6 0.8 1 

(b) e = 26.57° (100 x 100 subcells) 



0 0.5 1 1.5 


(c) 6= 14.04 (150 x 150 subcells) 



(d) 6 = 36.87° (150 x 150 subcells) 


Figure 4. Detailed volume discretizations of the four repeating unit cells in the rotated coordinate 
systems employed to accurately capture the geometric details within each unit cell. 


RUC Rotation Angle Subcell Discretization Fiber Volume Fraction 


tan 1 (j) 

= 0° 

100 x 100 

0.3468 

tan” 1 !!) 

= 14.04° 

150 x 150 

0.3476 

tan _1 (|) 

= 26.56° 

100 x 100 

0.3440 

tan _1 (|) 

= 36.87° 

150 x 150 

0.3511 

tan^ 1 (j) 

= 45° 

100 x 100 

0.3504 


Table 2. Geometric and microstructural details of the investigated repeating unit cells. 
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the use of square subcells to approximate the fiber shape. The actual subcell dimensions in each 
unit cell were dictated by the need to capture the circular fiber shape with sufficient detail given 
the required number of fibers within each repeating unit cell and the targeted fiber volume fraction. 
In the case of the repeating unit cell in the principal material coordinate system with a single fiber 
in the center (not shown), the repeating unit cell was discretized into 100 x 100 subcells in the 
— %3 plane. The same number was used for the repeating unit cells rotated by 26.57° and 45.0° 
about the x\ axis, while 150 x 150 subcells were employed for the unit cells rotated by 14.04° and 
36.87°. 

The repeating unit cell in the principal material coordinate system ( 6 = 0°, see the repeating 
unit cell with a single fiber in the center in Fig. 3) produces the homogenized elastic stiffness matrix 
C* of the form 
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(74) 


where C\ 2 = C* 3 , C 2 2 = C33, C g 5 = due to the cubic symmetry, but C44 ^ l,(C 22 — C 2 3) due 
to the absence of isotropy in the X 2 — £3 plane. These principal stiffness matrix elements can be 
used in the transformation equations 

C*(6) = TiC^T^ 1 (75) 

to generate the corresponding homogenized stiffness matrix elements of the effective medium in the 
rotated coordinate system independently of the micromechanics-based solution for the homogenized 
stiffness matrix of a repeating unit cell in the rotated coordinate system. The transformation 
matrices Ti and T2 for rotation by the angle 6 about the x\ axis are, 
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where m = cos# and n = sin#. They relate stress and strain quantities in the principal coordinate 
system, <x and e. to the corresponding quantities in the rotated (primed) coordinate system, a 
and e 
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a = Ti<t, 


£ — T2£ 


and are used to derive Eq. (75) from Hooke’s law in the principal material coordinate system. 
Under the above transformation, the homogenized stiffness matrix C *(0) acquires the following 
form in the rotated coordinate system 
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(76) 


where C{ 2 = C{ 3 , C| 2 = C&, C* 55 = C* 66 , and C 2 * 4 = -C 3 * 4 . 

The knowledge of the effective stiffness matrix elements in the rotated coordinate system makes 
it possible to generate the corresponding transformed compliance matrix elements from the inverse 
relationship 


S*(0) = [C*(0)]” 1 


(77) 


The transformed compliance matrix S*(0) has the same structural form as the stiffness matrix 
C *(9), with the individual elements C*j(9) replaced by the corresponding elements S*j(9). These 
elements, in turn, can be used to determine the effective engineering properties in the rotated 


coordinate system 

as follows 
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(78) 


where E* t (9) are the three Young’s moduli, u*j(9) (i ^ j ) are the major Poisson’s ratios, G*j(0) 
i ^ j ) are the three shear moduli, and 77 * 23 ; ^23 i are Lekhnitskii’s coefficients of mutual influence 
of the 1st and 2nd kind, respectively, Lekhnitskii [30]. These coefficients provide a measure of the 
extent of anisotropy in the X 2 — £3 plane introduced by the rotation angle 9 about the fiber axis. 
The coefficients of the 1st kind represent ratios of transverse normal to transverse shear strains 
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due to transverse shear loading only. Similarly, the coefficients of the 2nd kind represent ratios of 
transverse shear to transverse normal strains due to transverse normal loading only. 

In the following sections, we first determine the stiffness matrix elements of the two unidirec- 
tional composite systems in the principal material coordinate system by performing the high-fidelity 
micromechanical analysis on the simple repeating unit cell with the single fiber in the center, Fig. 3. 
We then use the values of these stiffness matrix elements, which have the form given by Eq. (74), in 
the transformation equations, Eqs. (75), to determine the homogenized stiffness matrix elements in 
the coordinate system rotated by the angle 9 about the fiber axis in the range 0° < 6 < 45°. These 
are then compared with the corresponding stiffness matrix elements obtained independently from 
the micromechanical analyses of the four remaining repeating unit cells in the rotated coordinate 
systems shown in Figs. 3 and 4. The determination of the engineering ’’constants” is then obtained 
by first determining the homogenized elastic compliance elements from Eq. (77), and then using 
the results in the definitions for engineering properties given by Eqs. (78). 

In order to further demonstrate the need for the development of the high-fidelity micromechanics 
model, we compare the results generated in the manner described above with the corresponding 
results obtained from the re-formulated version of the original GMC model. The reformulated 
version is required due to the large number of subcells used in the constructed repeating unit cells. 
In addition, we compare the local stress fields under normal and shear loading in the different 
repeating unit cells predicted by the high-fidelity model with the original model. 

4.1 Glass/Epoxy Unidirectional Composite 

Figure 5 illustrates the dependence of the elastic stiffness matrix elements C*j on the rotation angle 
6 for the glass/epoxy system with the Young’s moduli ratio Ef/E m = 20. The predictions generated 
by the high-fidelity model for the non-zero elastic stiffness matrix elements of the repeating unit 
cells shown in Fig. 4 have been normalized by the corresponding stiffness matrix elements C*j 
obtained from the repeating unit cell in the principal material coordinate system (with 6 = 0°), 
with the exception of the 4 and 4 elements which have been normalized by C * 4 . The predictions 
of the original GMC model, which have also been normalized by the corresponding elements C*j 
obtained from the high-fidelity model, are included in the figure. The actual values are given in 
the Appendix. As observed in the transformation equation results, while C 42 , and Cg 6 are 
insensitive to the rotation angle 9, the elements C| 2 , C| 3 , and C 44 exhibit substantial dependence, 
with and C 44 increasing as much as 15% and 30%, respectively, as 9 increases from 0° to 45°, 
and C 22 decreasing by more than 5%. The dependence of the coupling elements C 24 and C34 on 
the rotation angle is also dramatic, albeit the magnitudes are much smaller relative to C 4l . 

The high-fidelity model accurately captures the correct 9 dependence of the effective stiffness 
matrix elements, noting that the vertical scale ranges in the individual figures have been deliberately 
chosen so as to highlight the differences between the micromechanics-based predictions and the 
transformation equations. The slight departures from the transformation equation predictions are 
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O HFGMC 
□ GMC 

Transformation 


0.95 • 
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0 10 20 30 40 50 

( a ) C* 1 (0°) / C*, 




Figure 5. Normalized effective stiffness matrix elements C*j of the glass/epoxy unidirectional com- 
posite as a function of the rotation angle 9 about the fiber axis. Comparison of the original and 
the High-Fidelity GMC predictions with the transformation equations. 
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Figure 5. Continued. 


most likely due to the small variations in the repeating unit cell fiber volume fractions and the 
somewhat different circular fiber approximations by the employed square subcell discretizations. It 
is remarkable, however, that the C| 4 and (7| 4 behavior is captured so well in view of the relatively 
small magnitudes of these normal-shear coupling moduli. In contrast, the original GMC model 
is incapable of capturing the correct 9 dependence of the elastic stiffness matrix elements with 
the exception of C 41 which is theoretically insensitive to the rotation angle. In particular, the 
dependence of (7| 2 is predicted to rapidly decrease with increasing 9 and then increase in the 
considered range, which contrasts with the monotonically decreasing behavior exhibited by both 
the transformation equations and the high-fidelity model. Further, the actual magnitudes of this 
stiffness matrix element at the different off-axis angles are substantially lower relative to the high- 
fidelity model results, despite similar values for the repeating unit cell in the principal coordinate 
system ( 9 = 0°). Similar anomalous behavior is observed in the remaining elements of the stiffness 
matrix which is further discussed below in the context of the engineering moduli. Of particular 
importance is the original GMC model’s failure to correctly capture the anisotropic behavior due to 
the rotation angle manifested by the non-zero values of C 24 and (7| 4 . The original model predicts 
that these elements are identically zero for all rotation angles, effectively suggesting that orthotropic 
behavior is preserved in the rotated coordinate systems. 

The engineering moduli dependence on the rotation angle 9, calculated using Eqs. (77) and 
(78), for these repeating unit cells is illustrated in Fig. 6. As observed, the correlation between the 
high-fidelity model and transformation equation predictions is remarkable for the three Young’s 
and shear moduli, and the three major Poisson’s ratios. As expected, the axial Young’s modulus 
E 41 (9) remains nearly constant for the differently-oriented repeating unit cells, as does the major 
Poisson’s ratio ^12 (9) and the out-of-plane shear modulus G* l2 {9). The in-plane moduli, on the 
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Figure 6. Normalized engineering moduli of the glass/epoxy unidirectional composite as a function 
of the rotation angle 6 about the fiber axis. Comparison of the original and the High-Fidelity GMC 
predictions with the transformation equations. 
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Figure 6. Continued. 

other hand, exhibit substantial 6 dependence. While the transverse modulus -E| 2 ($) decreases with 
increasing rotation angle, the transverse Poisson’s ratio u 23 (0) and the transverse shear modulus 
G 23 (#) increase. For the rotation angle 9 = 45°, the decrease in the transverse Young’s modulus is 
more than 15% of the principal material coordinate system value. For the same rotation angle, the 
increase in the transverse shear modulus and Poisson’s ratio is more dramatic, with increases of 
approximately 30% and 25%, respectively. Equally remarkable is the correlation for the coefficients 
of mutual influence ? 7 223 , V 3 , 23 i ^ 23 , 2 > an d ^23 3 which couple the normal and shear responses. These 
coefficients give a measure of the extent of anisotropy in the X 2 — £3 plane caused by the absence 
of material planes of symmetry. 

In contrast, the predictions of the original GMC model are quite poor, with the exception of 
the Young’s modulus EJ' 1 (6 ) ), as suggested by the results in Fig. 5. In particular, the transverse 
modulus E| 2 (0) in the low off-axis angle range differs from the transformation and high-fidelity 
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model predictions by over 20 % which is greater than the difference in the corresponding stiffness 
matrix element C 2 2 in view of the pronounced differences in the transverse Poisson’s ratio v 23 W 
in the low off-axis angle. The differences in the out-of- plane shear modulus G\ 2 (9) are also on the 
order of 20 % in the entire off-axis range excluding 0 = 0 °. An even greater difference is obtained 
for the transverse shear modulus G 23 (9) for rotation angles greater than 30 ° due to the original 
model’s failure to predict any variation for all rotation angles. The differences between the high- 
fidelity and the original models is rooted in the absence of shear coupling due to the use of a 
linear displacement field in the original model, which results in constant strain and stress fields 
in the individual subcells. As discussed in the Introduction, the application of traction continuity 
conditions in an average sense between individual subcells in each row and column of subcells 
renders the corresponding traction components constant and their magnitudes are dictated by the 
most compliant subcell. In the case of the transverse shear stress this effect is further enhanced 
and results in a uniform shear stress throughout the entire repeating unit cell. The comparison of 
microscale stress fields discussed below illustrates this point more clearly. 

To illustrate the influence of shear coupling, microscale stress distributions predicted by the 
high-fidelity and original models are compared for the repeating unit cell rotated by 8 = 26 . 57 ° 
about the fiber axis and subjected to axial and transverse normal loading, and axial and transverse 
shear loading. Similar results have been observed for other rotation angles. Figure 7 compares 
the microscale cjh, 022 and <723 stress fields in this repeating unit cell for uniaxial loading by the 
average axial normal strain E n = 0 . 1 %, with the remaining faces of the repeating unit cell traction- 
free in the average sense. As observed, the axial < 7 n stress distributions predicted by both models 
are nearly identical, thereby producing essentially the same values of the axial Young’s modulus. 
This is not the case for the 0-22 and (723 stress distributions. In particular, the (J22 distribution 
predicted by original GMC model consists of parallel strips along the y 2 direction of the repeating 
unit cell, with the strips alternating in sign such that the average stress 022 remains zero due to 
equilibrium considerations. The parallel strip pattern arises due to the fact that the normal stress 
(722 is a traction component in each row of subcells along the y2 direction, which must remain 
constant in that row due to the linear displacement field approximation in each subcell and the 
imposition of traction continuity in a surface-averaged sense across subcell interfaces. The actual 
magnitude of (722 in each row is dictated by the relative content of the individual phases in that row, 
which represents the actual fiber arrangement only in a very approximate manner. This is a direct 
consequence of the lack of shear coupling and clearly reveals GMC to be a spring-like model with 
coupling only among the normal stress (strain) components. The absence of shear coupling results 
in the vanishing of the transverse shear stress (723 predicted by GMC in view of the vanishing of 
the average transverse shear stress on the external faces of each boundary subcell. In contrast, the 
high-fidelity model predictions of the microscale <722 and <723 stress distributions are consistent with 
the actual fiber arrangement. In particular, the distributions of both stress components within the 
individual fibers exhibit small departures from uniform distributions (as would be expected for the 
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Figure 7. Comparison of the an, 022 and <723 stress fields within the repeating unit cell of the 
glass/epoxy composite rotated by 26.57° about the fiber axis and subjected to the average axial 
normal strain en = 0.1%. High-Fidelity (left column) vs original GMC (right column) predictions. 
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dilute case) , with highly nonuniform distributions in the matrix phase which follows the geometric 
fiber pattern. Substantially higher magnitudes of these stress components relative to the GMC 
predictions, characterized by localized stress concentrations along the fiber/matrix interfaces, are 
observed due to the fiber/matrix elastic moduli mismatch and geometric fiber arrangement. It is 
remarkable that the transverse normal and shear stress components, induced by the application 
of just the axial stress, are reproduced with such good fidelity despite their substantially smaller 
magnitudes relative to the axial stress an. 

Figure 8 compares the microscale a 12 and <713 stress fields in this repeating unit cell for loading 
by the average axial shear strain £12 = 0.1%, with the remaining faces of the repeating unit 
cell traction- free in the average sense. In the case of the <712 stress component, the high-fidelity 
model predicts small departures from uniform stress distributions within the individual fibers, 
with magnitudes substantially greater than in the surrounding matrix due to the large fiber/matrix 
moduli mismatch. High stress concentrations are evident along certain segments of the fiber/matrix 
interfaces which may or may not be due to the step-wise discretization of the circular interface. 
This aspect requires further investigation which is beyond the scope of the present study (see 
Bednarcyk et al. [ 25 ] for related discussion about the mesh sensitivity of HFGMC). The ai3 stress 
distributions within the individual fibers also exhibit small departures from uniform distributions, 
although their magnitudes are nearly zero. In the matrix phase, however, this stress component is 
not insignificant relative to <712- In contrast, the ai2 distribution predicted by the original GMC 
model consists of a parallel strip pattern along the 1/2 axis, as was the case with the a22 distribution 
for loading by the axial normal stress. Significantly lower magnitudes, with little variation from 
strip to strip, are observed in the <712 distribution due to the weakest link effect caused by the 
application of the traction continuity conditions in an average sense in the presence of constant 
subcell stress held. This produces significantly lower axial shear modulus as already observed in 
Fig. 6 . Further, while the high-fidelity model predicts the presence of non-zero <713, this stress 
component is zero according to the original model. 

Figure 9 compares the microscale a 22, 033 and a 23 stress fields in this repeating unit cell 
for loading by the average transverse normal strain £22 = 0 . 1 %, with the remaining faces of the 
repeating unit cell traction-free in the average sense. The 022 stress distribution predicted by the 
high-fidelity model is qualitatively similar to the <712 distribution of the preceding example, while 
the a33 stress distribution is similar to a 13. Significant a 33 stress concentrations are present in the 
matrix phase at the fiber/matrix interfaces at points along the fiber diametral planes lined up with 
the load axis. Significant a23 stress magnitudes are also evident in oval regions surrounded by four 
fibers aligned with the rotated fiber rows. In contrast, and as expected from the previous results, 
the normal <722 and (J33 stress distributions generated by the original GMC model exhibit parallel 
strip patterns along the ij2 and y% directions, respectively. Little variation in the respective stress 
magnitudes is observed in the adjacent strips and the low magnitudes relative to the high-fidelity 
results produce a low value of the average normal stress 022 , thereby resulting in a low value of the 
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Figure 8. Comparison of the a\2 and a 13 stress fields within the repeating unit cell of the glass/epoxy 
composite rotated by 26 . 57 ° about the fiber axis and subjected to the average axial shear strain 
e\2 = 0 . 1 %. High-Fidelity (left column) vs original GMC (right column) predictions. 


transverse Young’s modulus for this rotation angle observed in Fig. 6 . Further, the transverse 
shear stress (723 is identically zero due to the absence of shear coupling. 

Finally, Figure 10 compares the microscale <723 and <722 stress fields in this repeating unit cell 
for loading by the average transverse shear strain £23 = 0 . 1 %, with the remaining faces of the 
repeating unit cell traction-free in the average sense. This loading case highlights the differences 
in the two stress distributions predicted by the original and the high-fidelity GMC models. The 
uniform (723 stress distribution throughout the entire repeating unit cell irrespective of location 
(i.e. , whether the particular point lies within the hard fiber or much softer matrix phase) predicted 
by the original GMC model is a direct consequence of the imposition of shear traction continuity 
in the surface-average sense across subcell interfaces in the \j2 and y 3 directions, given the linear 
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(e) ct 23 (HFGMC) (f) ct 23 (GMC) 


Figure 9 . Comparison of the U22, 033 and <723 stress fields within the repeating unit cell of the 
glass/epoxy composite rotated by 26 . 57 ° about the fiber axis and subjected to the average transverse 
normal strain £22 = 0 . 1 %. High-Fidelity (left column) vs original GMC (right column) predictions. 
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Figure 10. Comparison of the (723 and 022 stress fields within the repeating unit cell of the 
glass/epoxy composite rotated by 26.57° about the fiber axis and subjected to the average transverse 
shear strain £23 = 0.1%. High-Fidelity (left column) vs original GMC (right column) predictions. 



(c) a 22 (HFGMC) 




(a) a (HFGMC) 



displacement field approximation within each subcell. The low magnitude of this stress component 
relative to the high-fidelity prediction produces a substantially lower value of the transverse shear 
modulus observed in Fig. 6. The linear displacement field approximation also uncouples the 
shear and normal stress fields at the local level, thereby producing vanishing 022 stresses in the 
individual subcells. In contrast, the second-order displacement field approximation employed by 
the high-fidelity model is sufficient to correctly capture the stress transfer mechanism between the 
two phases, enabling the fibers to carry substantially higher shear stresses, thereby producing a 
higher average shear stress that the repeating unit cell can support for the same applied average 
shear strain. This results in the higher effective transverse shear modulus seen in Fig. 6. The 
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second-order displacement field approximation also couples the normal and shear stress fields at 
the local level, resulting in the observed 022 stress distribution which is nearly uniform in the 
individual fibers and highly nonuniform in the matrix phase. The matrix 0^2 stress nonuniformity 
is characterized by substantial stress concentrations at opposite locations around the fiber /matrix 
interface that are aligned with the diametral fiber planes coincident with the rotated fiber rows. 

4.2 Aluminum Weakened by Axially Oriented Cylindrical Porosities 

Figures 11 through 16 present the corresponding results for the aluminum matrix with substantially 
softer cylindrical inclusions, with the actual stiffness matrix and engineering moduli values given 
in the Appendix. In this case, the Young’s moduli ratio is Ef/E m = 0.01. For such a low ratio, 
the compliant cylindrical inclusions effectively behave as porosities. 

Figure 11 illustrates the dependence of the elastic stiffness matrix elements C* 3 on the rotation 
angle 8. The transformation equation predictions follow the trends presented in Fig. 5 for the 
glass/epoxy composite, but the variations are greater. In particular, the presence of porosities pro- 
duces a degradation in C | 2 greater than 15% as the rotation angle increases from 0° to 45°, which is 
roughly three times larger than the corresponding 5% decrease observed for the glass/epoxy system. 
The increases in the (7 | 3 and C 44 stiffness elements, which are more than 50% and nearly 100% 
respectively, are also substantially greater relative to the glass/epoxy system. The normal/shear 
coupling elements 4 and C34 also experience greater variations, with the maximum values be- 
ing approximately five-fold greater relative to the glass/epoxy composite. The high-fidelity GMC 
model predictions follow the transformation equation variations very closely, with the observed 
small departures most likely due to the small variations in the porosity content among the five 
different repeating unit cells containing different numbers of porosities. 

In contrast, the original GMC model is not capable of accurately predicting the magnitude and 
variation of any of the stiffness elements with the rotation angle. The determination of the individual 
stiffness matrix elements involves the application of one nonzero macroscopic strain component at a 
time which gives rise to tractions whose line of action traverses the porous or compliant phase even 
in the case of the axial strain. Under transverse loading, the stress cannot be transferred effectively 
around the compliant inclusions due to the absence of shear coupling, and the load supported by 
any row of subcells coincident with the applied load is governed by the most compliant subcell 
in that row. This produces substantial degradation in the values of the stiffness elements in the 
principal coordinate system, which is greatly compounded by the rotation angle, including the 
( 7*4 stiffness element which was so accurately predicted for the glass/epoxy composite. It is clear 
that for this porous material, the original GMC predictions for all stiffness matrix elements at any 
non-zero rotation angle are either significantly different or completely erroneous. It is only for the 
stiffness elements which directly relate the normal stress-strain pairs at the zero rotation angle, 
C 41 , G 22 , and the axial shear pairs, C 55 , Cg 6 , that the predictions are not dramatically different. In 
contrast, the stiffness element that relates the transverse shear stress-strain pair, C| 4 , is completely 
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Figure 11. Normalized effective stiffness matrix elements C*j of the aluminum matrix weakened by 
cylindrical porosities as a function of the rotation angle 8 about the porosity axis. Comparison of 
the original and the High-Fidelity GMC predictions with the transformation equations. 
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Figure 11. Continued. 


erroneous even at the zero rotation angle. This will be further examined in the Discussion section 
in the context of using an alternative method to determine the off-axis stiffness element behavior 
using the original GMC model. 

The greater variations in the elastic stiffness matrix elements with the rotation angle for this 
material system translate into correspondingly greater moduli variations shown in Fig. 12. As 
in the case of the glass/epoxy composite, the transformation equations predict variations only for 
the transverse moduli and the normal/shear coupling coefficients. In particular, the transverse 
modulus E 22 (9) decreases by nearly 40% at the rotation angle of 45° relative to its value in the 
principal material coordinate system, compared to just a little more than 15% for the glass/epoxy 
system. The increase in the transverse modulus G 23 (9) is almost 100%, which is the same as the 
increase in (7| 4 (0). The increase in the transverse Poisson’s ratio ^ 23 ( 0 ) is even more dramatic, 
being around 150%. These large variations with the rotation angle are captured very well by the 
high-fidelity GMC model, as expected from the stiffness matrix element results. In contrast, the 
original GMC model results are completely erroneous for almost all engineering moduli at non-zero 
rotation angles as also suggested by the preceding stiffness matrix element results. The exception 
is the axial Young’s modulus £'i 1 (0) and the related Poisson’s ratio v\ 2 which are predicted very 
accurately for all rotation angles. In the case of pure axial loading in the porosity direction, the 
matrix phase is continuous along this direction, and thus effective in supporting the entire axial 
load without the need for stress transfer through the shear-normal coupling mechanism, while the 
transverse contraction which affects v\ 2 occurs unconstrained due to porosity’s presence. 

The above results are explained by microscale stress distributions predicted by the high-fidelity 
and original models presented in Figs. 13 through 16. As in the case of the glass/epoxy composite, 
these distributions have been generated for the repeating unit cell rotated by 6 = 26.57° and 
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Figure 12. Normalized engineering moduli of the aluminum matrix weakened by cylindrical porosi- 
ties as a function of the rotation angle 9 about the porosity axis. Comparison of the original and 
the High-Fidelity GMC predictions with the transformation equations. 
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Figure 12 . Continued. 


subjected to axial and transverse normal loading, and transverse and axial shear loading. Figure 
13 compares the microscale an, a 22 and <723 stress fields in this repeating unit cell for loading by 
the average axial normal strain e\\ = 0 . 1 %, with the remaining faces of the repeating unit cell 
traction-free in the average sense. As in the preceding case, virtually no difference is observed 
in the an stress distributions predicted by both models. In contrast to the preceding case, no 
differences are observed in the 022 and a23 stress fields as well. These stress components are, in 
fact, very small due to the absence of constraint on the transverse deformation of the matrix under 
pure axial loading in the presence of porosities. This lack of constraint allows the matrix phase to 
contract freely without inducing transverse normal and shear stresses. This is one instance when 
the original GMC model correctly predicts the microscale stress fields which, in turn, produce the 
correct variation in the axial modulus ( 0 ). 
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Figure 13. Comparison of the an, 022 and a 23 stress fields within the repeating unit cell of porous 
aluminum matrix rotated by 26.57° about the porosity axis and subjected to the average axial 
normal strain £11 = 0.1%. High-Fidelity (left column) vs original GMC (right column) predictions. 
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Figure 14. Comparison of the a 12 and <713 stress fields within the repeating unit cell of porous 
aluminum matrix rotated by 26.57° about the porosity axis and subjected to the average axial 
shear strain £12 = 0.1%. High-Fidelity (left column) vs original GMC (right column) predictions. 

Figure 14 compares the microscale a 12 and <713 stress fields in this repeating unit cell for loading 
by the average axial shear strain £12 = 0.1%, with the remaining faces of the repeating unit cell 
traction- free in the average sense. While the compliant inclusions or porosities cannot support 
appreciable axial shear stresses, the matrix phase does support these stress components. This is 
correctly predicted by the high-fidelity model where appreciable levels of both stress components 
are carried by the matrix phase. The highly nonuniform stress distributions are a direct result of 
the shear coupling effect which is responsible for transferring the two stress components around 
porosities. This stress transfer mechanism is absent in the original GMC model, resulting in very 
low values of the axial shear stress <712 (noting that the axial shear modulus G* 2 ( 9 ) is not quite 
zero), and zero values of the axial shear stress <713 due to decoupling of the axial shear response in 
the yi — y2 and y\ — y 3 planes in the case of isotropic or orthotropic phases. 
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Figure 15 compares the microscale 022, 033 and 023 stress fields in this repeating unit cell 
for loading by the average transverse normal strain £22 = 0 . 1 %, with the remaining faces of the 
repeating unit cell traction-free in the average sense. The detrimental effect of the shear coupling’s 
absence in the original GMC model on these stress components is clearly observed, with the 022 
stress field characterized by essentially uniform and very low magnitudes, and thus a low transverse 
Young’s modulus £ l | 2 (^) seen i n Fig- 12 . The same holds true for the <733 stress field, while the 
transverse shear stress <723 vanishes completely. In contrast, the shear coupling effects necessary 
to internally support the applied transverse load are clearly evident in the high-fidelity GMC 
model predictions for the three stress fields. Highly nonuniform distributions are observed for the 
three stress components in the matrix phase, characterized by significant concentrations at the 
fiber/matrix interfaces at specific locations as in the preceding case. 

Similar trends in the stress distributions predicted by the two models are observed for loading 
by the average transverse shear strain £23 = 0 . 1 %, with the remaining faces of the repeating unit 
cell traction-free in the average sense. This is seen Fig. 16 which compares the microscale 023 and 
022 stress fields within the considered repeating unit cell. Interestingly, large concentrations of the 
normal stress <7 22 are observed at the fiber/matrix interface due to the transverse shear loading in 
the high-fidelity GMC model predictions. These are higher than the corresponding transverse shear 
stress concentrations. The original GMC model is incapable of capturing these nonuniform stress 
fields due to the absence of shear coupling, and the presence of porosities results in uniform and very 
low magnitudes of the transverse shear stress, producing a very low value of the transverse shear 
modulus G^iO) seen in Fig. 12 . Furthermore, the normal stress 022 which exhibits high stress 
concentrations observed in the high-fidelity result completely vanishes according to the original 
GMC model. 

5 DISCUSSION 

As illustrated in the foregoing, for a repeating unit cell without planes of material symmetry 
parallel to the fiber axis the microscale stress fields predicted by the original GMC model exhibit 
two characteristic patterns in the plane normal to the fiber axis. The normal <722 and 033 stress 
distributions are characterized by parallel strips whose magnitude and sign depends on the applied 
stress orientation and the fiber/matrix Young’s modulus ratio. For large ratios the pattern is 
pronounced while for very low ratios which mimic porosities the patterns are obscured by the very 
low normal stress magnitudes. Similar strip patterns are observed for repeating unit cells with a 
single fiber in the center, which possess two orthogonal planes of material symmetry, Bednarcyk 
et al. [ 25 ] and Pindera et al. [ 26 ]. In such cases, the strip patterns are wider and fewer due 
to the absence of overlapping fibers aligned with the load axis. In contrast, the shear 023 stress 
distributions are uniform within the rotated repeating unit cells as well as within unit cells in the 
principal material coordinate system. As discussed herein, and also elsewhere, these patterns are a 
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Figure 15 . Comparison of the <J22 j <733 and <723 stress fields within the repeating unit cell of porous 
aluminum matrix rotated by 26 . 57 ° about the porosity axis and subjected to the average transverse 
normal strain £22 = 0 . 1 %. High-Fidelity (left column) vs original GMC (right column) predictions. 
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Figure 16. Comparison of the a 23 and 022 stress fields within the repeating unit cell of porous 
aluminum matrix rotated by 26.57° about the porosity axis and subjected to the average transverse 
shear strain £23 = 0.1%. High-Fidelity (left column) vs original GMC (right column) predictions. 

direct result of the absence of shear coupling and lead to inaccurate microscale stress fields which 
cannot be used to predict local fiber or matrix failure caused by stress combinations other than 
the second stress deviator. As also shown for the porous aluminum case, these characteristic stress 
distributions produce highly inaccurate engineering moduli in the rotated coordinate system. 

In order to mitigate the negative impact of the shear coupling absence, an alternative manner 
of determining the repeating unit cell response, and thus the engineering moduli, in the rotated 
coordinate system based on the GMC model can be chosen. First, the applied normal or shear 
strain in the rotated coordinate system is transformed to the principal coordinate system in order 
to determine the response in this reference frame. Then the resulting strains and stresses are 
transformed back to the rotated coordinate system in order to determine the elastic moduli. In 
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fact, this is the basis for the transformation equations given by Eq. (75), and should yield the same 
result as direct application of strains in the rotated coordinate system if the micromechanics model 
is self consistent. This is clearly the case for the HFGMC model. In the case of GMC, however, 
the observed characteristic stress patterns indicate that this is going to be only partially successful 
in the presence of very compliant inclusions. In particular, transforming the applied macroscopic 
strain to the principal material coordinate system will not negate the fact that the transverse shear 
response, and thus the transverse shear modulus, will be incorrectly predicted in the presence 
of porosities. In particular, the transverse shear modulus will be zero rendering the resulting 
macroscopic shear stress zero. Therefore, this contribution will be absent when transforming the 
stresses to the rotated coordinate system, thereby producing a result that will differ from the 
transformation equations by an amount that depends on the rotation angle or the magnitude of the 
absent transverse shear stress in the principal material coordinate system. However, the dramatic 
differences between GMC predictions and transformation equations observed in Fig. 11 and 12 
for the porous aluminum case will be reduced because the contribution of the normal stresses in 
the principal material coordinate system will not be completely eliminated even in the presence of 
porosities. This method of calculating the engineering moduli in the rotated coordinate system is 
in fact equivalent to the use of transformation equations based on the moduli calculated by the 
GMC model in the principal material coordinate system. The different results obtained from the 
GMC-based calculations, which depend on whether the calculations are made in the rotated or 
unrotated coordinate system, point to a fundamental problem which limits this model’s range of 
applicability. This limitation has clearly been overcome by the high-fidelity version which can be 
used to accurately model the response of a wide range of periodic materials with or without planes 
of material symmetry in different coordinate systems. 

6 SUMMARY AND CONCLUSIONS 

The reformulation of the High-Fidelity Generalized Method of Cells, based on a simplified vol- 
ume discretization involving only subcells as the fundamental subvolumes together with the use 
of the local/global stiffness matrix approach, facilitates the analysis of repeating unit cells with 
complex microstructural details characteristic of realistic microstructures of multiphase materials. 
This is a direct result of the elimination of redundant continuity conditions present in the original 
formulation, which in turn produces substantial reduction in the size of the system of equations 
governing the response of the repeating unit cell. The reformulation also reveals the high-fidelity 
micromechanical analysis to be an approximate elasticity technique based on a direct volume av- 
eraging of the stress equilibrium equations and the imposition of the displacement and traction 
continuity conditions in a surface- averaged sense across the interfaces between adjacent subcells. 
This, in turn, simplifies the derivation of the volume-averaged equilibrium equations governing the 
individual subcell response as well as the derivation of the traction continuity conditions. 
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In the present investigation, the reformulation was employed to determine the elastic moduli of 
a square array of stiff fibers embedded in a substantially more compliant matrix, representative of a 
unidirectional glass/epoxy composite, under rotation about the fiber axis. A stiff matrix weakened 
by cylindrical porosities was also considered. The rotation about the fiber axis necessitates the 
analysis of repeating unit cells in the rotated coordinate system, representative of the same square 
array, which may contain many fibers, in contrast with the single fiber of the repeating unit cell in 
the principal material coordinate system. Such repeating unit cells typically do not possess planes 
of material symmetry, rendering them anisotropic in the rotated coordinate system. The elastic 
moduli of these unit cells can also be obtained from the standard transformation equations once 
the effective composite properties have been calculated in the principal material coordinate system. 
These transformation equations can therefore be employed to validate the predictive capability of a 
micromechanics model that admits periodic boundary conditions required in the absence of material 
planes of symmetry. The effective moduli predicted by the High-Fidelity Generalized Method of 
Cells based on the unit cells in the rotated coordinates systems have been shown herein to correlate 
extremely well with the transformation equations for both material systems considered. In contrast, 
the predictions of the original Generalized Method of Cells exhibited substantial departures from 
the transformation equations for the glass/epoxy system, which became unacceptably large for 
the porous aluminum. This is a direct result of the absence of shear coupling in the original 
method which produces erroneous results when normal/shear interaction dominates in the presence 
of porosities or inclusion phases that are substantially more compliant than the matrix phase. The 
high-fidelity version circumvents this problem, albeit at an increased computational cost which, 
however, is mitigated to a certain extent by the implemented reformulation. The microscale stress 
distributions generated by both models for different uniaxial loading situations provided additional 
insight supporting the predicted moduli results. 
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8 APPENDIX 


8.1 Local Stiffness Matrix Elements 

Explicit expressions for the non-zero elements of the local stiffness matrix for the ((3, 7 ) subcell, 
given in terms of the subcell’s geometric and mechanical properties, are listed below for loading by 
normal and shear tractions in the 1/2 — 1/3 plane. 
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CELL EFFECTIVE STIFFNESS 

2 . 9214E+04 4 . 2656E+03 

4 . 2656E+03 9.5857E+03 

4 . 2656E+03 4.1918E+03 

-4 . 6064E-11 -6 . 7785E-11 

0 . OOOOE+OO 0 . OOOOE+OO 
0 . OOOOE+OO 0. OOOOE+OO 

EFFECTIVE MODULI 

Ell = 2 . 6573E+04 

E22 = 7 . 5417E+03 

E33 = 7 . 5417E+03 

G23 = 2 . 0856E+03 

G13 = 2 . 5070E+03 

G12 = 2 . 5070E+03 

NU23 = 3 . 9820E-01 

NU13 = 3 . 0961E-01 

NU12 = 3 . 0961E-01 


4 . 2656E+03 -4.6827E-11 

4 . 1918E+03 -9 . 7567E-11 

9 . 5857E+03 -9.6209E-11 

-1 . 2299E-10 2 . 0856E+03 

0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 


0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
2 . 5070E+03 -4 . 4269E-11 

-4 . 3698E-11 2 . 5070E+03 
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ANGLE = 14.04 degrees 


CELL EFFECTIVE STIFFNESS 


2 . 9268E+04 

4 . 2738E+03 

4 . 2738E+03 

8 . 9032E-12 

0. OOOOE+OO 

0. OOOOE+OO 

4 . 2738E+03 

9 . 5079E+03 

4 . 3033E+03 

-2 . 5483E+02 

0. OOOOE+OO 

0. OOOOE+OO 

4 . 2738E+03 

4 . 3033E+03 

9 . 5079E+03 

2 . 5483E+02 

0. OOOOE+OO 

0. OOOOE+OO 

1 . 0230E-11 

-2 . 5483E+02 

2 . 5483E+02 

2 . 2148E+03 

0. OOOOE+OO 

0. OOOOE+OO 

0 . OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

2 . 5218E+03 

1 . 9721E-11 

0 . OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

1 . 9821E-11 

2 . 5218E+03 


EFFECTIVE MODULI 

Ell = 2 . 6623E+04 
E22 = 7 . 3013E+03 
E33 = 7 . 3013E+03 
G23 = 2 . 1899E+03 
G13 = 2 . 5218E+03 
G12 = 2 . 5218E+03 
NU23 = 4 . 1883E-01 
NU13 = 3 . 0944E-01 
NU12 = 3 . 0944E-01 


ANGLE = 26.57 degrees 


CELL EFFECTIVE STIFFNESS 


2 . 9023E+04 

4 . 2544E+03 

4 . 2544E+03 

5 . 1096E-12 

0. OOOOE+OO 

0. OOOOE+OO 

4 . 2544E+03 

9 . 1527E+03 

4 . 5782E+03 

-2 . 8695E+02 

0. OOOOE+OO 

0. OOOOE+OO 

4 . 2544E+03 

4 . 5782E+03 

9 . 1527E+03 

2 . 8695E+02 

0. OOOOE+OO 

0. OOOOE+OO 

7 . 6872E-12 

-2 . 8695E+02 

2 . 8695E+02 

2.4771E+03 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

2 . 5013E+03 

5 . 9124E-11 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

5 . 8265E-11 

2 . 5013E+03 


EFFECTIVE MODULI 

Ell = 2 . 6387E+04 
E22 = 6 . 6240E+03 
E33 = 6 . 6240E+03 
G23 = 2 . 4411E+03 
G13 = 2 . 5013E+03 
G12 = 2 . 5013E+03 
NU23 = 4 . 6938E-01 
NU13 = 3 . 0984E-01 
NU12 = 3 . 0984E-01 


ANGLE = 36.87 degrees 


CELL EFFECTIVE STIFFNESS 


2 . 9514E+04 

4 . 2947E+03 

4 . 2947E+03 

-5 . 4990E-12 

0. OOOOE+OO 

0. OOOOE+OO 

4 . 2947E+03 

9 . 1142E+03 

4 . 7835E+03 

-1 . 7227E+02 

0. OOOOE+OO 

0. OOOOE+OO 

4 . 2947E+03 

4 . 7835E+03 

9 . 1142E+03 

1 . 7227E+02 

0. OOOOE+OO 

0. OOOOE+OO 

-5 . 5649E-12 

-1 . 7227E+02 

1 . 7227E+02 

2 . 6792E+03 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

2 . 5480E+03 

6 . 7261E-02 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

7 . 7051E-12 

2 . 5480E+03 
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EFFECTIVE MODULI 


Ell = 

2 . 6860E+04 

E22 = 

6 . 4275E+03 

E33 = 

6 . 4275E+03 

G23 = 

2 . 6655E+03 

G13 = 

2 . 5480E+03 

G12 = 

2 . 5480E+03 

NU23 = 

4 . 9181E-01 

NU13 = 

3 . 0902E-01 

NU12 = 

3 . 0902E-01 


ANGLE = 45 degrees 


CELL EFFECTIVE STIFFNESS 


2 . 9462E+04 

4 . 2848E+03 

4 . 2848E+03 

-3 . 1417E-11 

0 . 0000E+00 

0 . 0000E+00 

4 . 2848E+03 

9 . 0188E+03 

4 . 8382E+03 

-3 . 5367E-11 

0 . 0000E+00 

0 . 0000E+00 

4 . 2848E+03 

4 . 8382E+03 

9 . 0188E+03 

-9 . 4672E-11 

0 . 0000E+00 

0 . 0000E+00 

-2 . 6528E-11 

-1 . 9263E-11 

-9 . 1672E-11 

2 . 7397E+03 

0 . 0000E+00 

0 . 0000E+00 

0 . 0000E+00 

0 . 0000E+00 

0 . 0000E+00 

0 . 0000E+00 

2 . 5268E+03 

3 . 4886E-11 

0 . 0000E+00 

0 . 0000E+00 

0 . 0000E+00 

0 . 0000E+00 

3.4856E-11 

2 . 5268E+03 


EFFECTIVE MODULI 

Ell = 2 . 6812E+04 
E22 = 6 . 2794E+03 
E33 = 6 . 2794E+03 
G23 = 2 . 7397E+03 
G13 = 2 . 5268E+03 
G12 = 2 . 5268E+03 
NU23 = 5 . 0206E-01 
NU13 = 3 . 0922E-01 
NU12 = 3 . 0922E-01 


8.2.2 Glass/Epoxy Composite: GMC Predictions 


ANGLE = 0 degrees 


CELL EFFECTIVE STIFFNESS 


0 . 2918E+05 
0 . 4186E+04 
0 . 4186E+04 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


0 . 4186E+04 
0 . 9326E+04 
0 . 4124E+04 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


0 . 4186E+04 
0 . 4124E+04 
0 . 9326E+04 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 1937E+04 
0 . 0000E+00 
0 . 0000E+00 


0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 2266E+04 
0 . 0000E+00 


0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 2266E+04 


EFFECTIVE MODULI 


El 1S= 0 . 2657E+05 
N12S= 0 . 3113E+00 
E22S= 0 . 7303E+04 
N23S= 0 . 4037E+00 
E33S= 0 . 7303E+04 
G23S= 0 . 1937E+04 
G13S= 0 . 2266E+04 
G12S= 0 . 2266E+04 
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ANGLE = 14.04 degrees 


CELL EFFECTIVE STIFFNESS 


0 . 2914E+05 
0 . 4015E+04 
0 . 4015E+04 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


0 . 4015E+04 
0 . 8318E+04 
0 . 4422E+04 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


0 . 4015E+04 
0 . 4422E+04 
0 . 8318E+04 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 


0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 1939E+04 
0 . 0000E+00 
0 . 0000E+00 


0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 0000E+00 
0 . 1945E+04 
0 . 0000E+00 


EFFECTIVE MODULI 


E11S= 0 . 2661E+05 
N12S= 0 . 3151E+00 
E22S= 0 . 5838E+04 
N23S= 0 . 4982E+00 
E33S= 0 . 5838E+04 
G23S= 0 . 1939E+04 
G13S= 0 . 1945E+04 
G12S= 0 . 1945E+04 


ANGLE = 26.57 degrees 


CELL EFFECTIVE STIFFNESS 


0 . 2890E+05 
0 . 4005E+04 
0 . 4005E+04 
0 . OOOOE+OO 
0 . OOOOE+OO 
0. OOOOE+OO 


0 . 4005E+04 
0 . 8298E+04 
0 . 4400E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 4005E+04 
0 . 4400E+04 
0 . 8298E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1929E+04 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1942E+04 
0. OOOOE+OO 


EFFECTIVE MODULI 


El 1S= 0 . 2638E+05 
N12S= 0 . 3154E+00 
E22S= 0 . 5834E+04 
N23S= 0 . 4966E+00 
E33S= 0 . 5834E+04 
G23S= 0 . 1929E+04 
G13S= 0 . 1942E+04 
G12S= 0 . 1942E+04 


ANGLE = 36.87 degrees 


CELL EFFECTIVE STIFFNESS 


0 . 2938E+05 
0 . 4026E+04 
0 . 4026E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 4026E+04 
0 . 8342E+04 
0 . 4443E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 4026E+04 
0 . 4443E+04 
0 . 8342E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1949E+04 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1949E+04 
0. OOOOE+OO 


EFFECTIVE MODULI 
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0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1945E+04 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1942E+04 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1949E+04 



El 1S= 0 . 2685E+05 
N12S= 0 . 3149E+00 
E22S= 0 . 5846E+04 
N23S= 0 . 4996E+00 
E33S= 0 . 5846E+04 
G23S= 0 . 1949E+04 
G13S= 0 . 1949E+04 
G12S= 0 . 1949E+04 


ANGLE = 45 degrees 


CELL EFFECTIVE STIFFNESS 


0 . 2936E+05 
0 . 4064E+04 
0 . 4064E+04 
0 . OOOOE+OO 
0 . OOOOE+OO 
0. OOOOE+OO 


0 . 4064E+04 
0 . 8511E+04 
0 . 4430E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 4064E+04 
0 . 4430E+04 
0 . 8511E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1947E+04 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 2008E+04 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 2008E+04 


EFFECTIVE MODULI 


El 1S= 0 . 2680E+05 
N12S= 0 . 3140E+00 
E22S= 0 . 6066E+04 
N23S= 0 . 4866E+00 
E33S= 0 . 6066E+04 
G23S= 0 . 1947E+04 
G13S= 0 . 2008E+04 
G12S= 0 . 2008E+04 


8.2.3 Aluminum with Cylindrical Porosities: HFGMC Predictions 

ANGLE = 0 degrees 


CELL EFFECTIVE STIFFNESS 


5 . 7040E+04 
1 . 6777E+04 
1 . 6777E+04 
-1 . 0418E-11 
0. OOOOE+OO 
0. OOOOE+OO 


1 . 67T7E+04 
3 . 8902E+04 
1 . 1938E+04 
-1 . 6283E-11 
0. OOOOE+OO 
0. OOOOE+OO 


1 . 6777E+04 
1 . 1938E+04 
3 . 8902E+04 
-2 . 1946E-11 
0. OOOOE+OO 
0. OOOOE+OO 


-2 . 4291E-11 
-1 . 9329E-11 
-5 . 3848E-11 
6 . 9500E+03 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
1 . 2819E+04 
7 . 5422E-11 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
7 . 5322E-11 
1 . 2819E+04 


EFFECTIVE MODULI 

Ell = 4 . 5967E+04 
E22 = 3 . 2523E+04 
E33 = 3 . 2523E+04 
G23 = 6 . 9500E+03 
G13 = 1 . 2819E+04 
G12 = 1 . 2819E+04 
NU23 = 2 . 0618E-01 
NU13 = 3 . 3000E-01 
NU12 = 3 . 3000E-01 
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ANGLE = 14.04 degrees 


CELL EFFECTIVE STIFFNESS 

5 . 6619E+04 1 . 6219E+04 
1 . 6219E+04 3 . 6558E+04 
1 . 6219E+04 1 . 2591E+04 
1 . 7596E-12 -2 . 6832E+03 
0 . 0000E+00 0 . 0000E+00 
0 . 0000E+00 0 . 0000E+00 

EFFECTIVE MODULI 

Ell = 4 . 5914E+04 
E22 = 2 . 8437E+04 
E33 = 2 . 8437E+04 
G23 = 7 . 1315E+03 
G13 = 1 . 2625E+04 
G12 = 1 . 2625E+04 
NU23 = 2 . 8651E-01 
NU13 = 3 . 3000E-01 
NU12 = 3 . 3000E-01 


ANGLE = 26.57 degrees 


CELL EFFECTIVE STIFFNESS 

5 . 7057E+04 1.6510E+04 

1 . 6510E+04 3 . 4149E+04 

1 . 6510E+04 1 . 5880E+04 

-2 . 1006E-11 -3 . 1257E+03 

0 . OOOOE+OO 0 . OOOOE+OO 
0 . OOOOE+OO 0. OOOOE+OO 

EFFECTIVE MODULI 

Ell = 4 . 6161E+04 

E22 = 2 . 3445E+04 

E33 = 2 . 3445E+04 

G23 = 9 . 9834E+03 

G13 = 1 . 2762E+04 

G12 = 1 . 2762E+04 

NU23 = 4 . 2077E-01 

NU13 = 3 . 3000E-01 

NU12 = 3 . 3000E-01 


ANGLE = 36.87 degrees 


CELL EFFECTIVE STIFFNESS 

5 . 6009E+04 1 . 5668E+04 

1 . 5668E+04 3 . 0781E+04 

1 . 5668E+04 1 . 6698E+04 

-4 . 6593E-12 -1.8320E+03 

0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 

EFFECTIVE MODULI 


1 . 6219E+04 3 . 4135E-12 
1 . 2591E+04 -2 . 6832E+03 
3 . 6558E+04 2.6832E+03 
2 . 6832E+03 7.7323E+03 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 


1 . 6510E+04 -1 . 6233E-11 
1 . 5880E+04 -3 . 1257E+03 
3 . 4149E+04 3 . 1257E+03 
3 . 1257E+03 1 . 1053E+04 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 


1 . 5668E+04 -6 . 1709E-12 
1 . 6698E+04 -1 . 8320E+03 
3 . 0781E+04 1 . 8320E+03 
1 . 8320E+03 1 . 2302E+04 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 


0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
1 . 2625E+04 -1 . 1416E-11 

-1.1035E-11 1 . 2625E+04 


0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
1 . 2762E+04 -7 . 1498E-13 

-8 . 7295E-13 1.2762E+04 


0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
0. OOOOE+OO 0. OOOOE+OO 
1 . 2408E+04 -1 . 3313E-12 

-1 . 0340E-12 1 . 2408E+04 
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Ell = 

4 . 5668E+04 

E22 = 

2 . 0060E+04 

E33 = 

2 . 0060E+04 

G23 = 

1 . 1825E+04 

G13 = 

1 . 2408E+04 

G12 = 

1 . 2408E+04 

NU23 = 

4 . 8183E-01 

NU13 = 

3 . 3000E-01 

NU12 = 

3 . 3000E-01 


ANGLE = 45 degrees 


CELL EFFECTIVE STIFFNESS 


5 . 6602E+04 
1 . 6492E+04 
1 . 6492E+04 
-8 . 2076E-12 
0 . OOOOE+OO 
0 . OOOOE+OO 


1 . 6492E+04 
3 . 1628E+04 
1 . 8348E+04 
4 . 6414E-12 
0. OOOOE+OO 
0. OOOOE+OO 


1 . 6492E+04 
1 . 8348E+04 
3 . 1628E+04 
-2 . 5772E-11 
0. OOOOE+OO 
0. OOOOE+OO 


-6 . 5033E-12 
2 . 3342E-11 
-4 . 2766E-11 
1 . 3374E+04 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
1 . 2681E+04 
4 . 1566E-12 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
3 . 5532E-12 
1 . 2681E+04 


EFFECTIVE MODULI 


Ell = 4.5717E+04 
E22 = 1 . 9985E+04 
E33 = 1 . 9985E+04 
G23 = 1 . 3374E+04 
G13 = 1 . 2681E+04 
G12 = 1 . 2681E+04 
NU23 = 5 . 0491E-01 
NU13 = 3 . 3000E-01 
NU12 = 3 . 3000E-01 


8.2.4 Aluminum with Cylindrical Porosities: GMC Predictions 


ANGLE = 0 degrees 


CELL EFFECTIVE STIFFNESS 


0 . 5359E+05 
0 . 1156E+05 
0 . 1156E+05 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 1156E+05 
0 . 2925E+05 
0 . 5778E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 1156E+05 
0 . 5778E+04 
0 . 2925E+05 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 7448E+03 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 9319E+04 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 9319E+04 


EFFECTIVE MODULI 


El 1S= 0 . 4597E+05 
N12S= 0 . 3300E+00 
E22S= 0 . 2635E+05 
N23S= 0 . 1228E+00 
E33S= 0 . 2635E+05 
G23S= 0 . 7448E+03 
G13S= 0 . 9319E+04 
G12S= 0 . 9319E+04 
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ANGLE = 14.04 degrees 

CELL EFFECTIVE STIFFNESS 

0 . 4687E+05 0.1452E+04 0.1452E+04 

0 . 1452E+04 0 . 2958E+04 0.1442E+04 

0 . 1452E+04 0 . 1442E+04 0.2958E+04 

0 . 0000E+00 0 . 0000E+00 0.0000E+00 

0 . 0000E+00 0 . 0000E+00 0.0000E+00 

0 . 0000E+00 0 . 0000E+00 0.0000E+00 

EFFECTIVE MODULI 

E11S= 0 . 4591E+05 

N12S= 0 . 3300E+00 
E22S= 0 . 2243E+04 
N23S= 0 . 4796E+00 
E33S= 0 . 2243E+04 
G23S= 0 . 7432E+03 
G13S= 0 . 7528E+03 
G12S= 0 . 7528E+03 


ANGLE = 26.57 degrees 

CELL EFFECTIVE STIFFNESS 

0.4713E+05 0 . 1475E+04 0.1475E+04 

0 . 1475E+04 0 . 3016E+04 0.1454E+04 

0 . 1475E+04 0 . 1454E+04 0.3016E+04 

0 . OOOOE+OO 0 . OOOOE+OO O.OOOOE+OO 
0 . OOOOE+OO O.OOOOE+OO 0.0000E+00 
0.0000E+00 O.OOOOE+OO O.OOOOE+OO 

EFFECTIVE MODULI 

E11S= 0 . 4616E+05 
N12S= 0 . 3300E+00 
E22S= 0 . 2303E+04 
N23S= 0 . 4740E+00 
E33S= 0 . 2303E+04 
G23S= 0 . 7507E+03 
G13S= 0 . 7699E+03 
G12S= 0 . 7699E+03 


ANGLE = 36.87 degrees 

CELL EFFECTIVE STIFFNESS 

0 . 4661E+05 0 . 1429E+04 0.1429E+04 

0 . 1429E+04 0 . 2901E+04 0.1429E+04 

0 . 1429E+04 0 . 1429E+04 0.2901E+04 

O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 

EFFECTIVE MODULI 

E11S= 0 . 4567E+05 


O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
0 . 7432E+03 O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO 0 . 7528E+03 O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO 0.7528E+03 


O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
0 . 7507E+03 O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO 0 . 7699E+03 O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO 0.7699E+03 


O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 
0 . 7359E+03 O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO 0 . 7361E+03 O.OOOOE+OO 
O.OOOOE+OO O.OOOOE+OO 0.7361E+03 
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N12S= 0 . 3300E+00 
E22S= 0 . 2186E+04 
N23S= 0 . 4847E+00 
E33S= 0 . 2186E+04 
G23S= 0 . 7359E+03 
G13S= 0 . 7361E+03 
G12S= 0 . 7361E+03 


ANGLE = 45 degrees 
CELL EFFECTIVE STIFFNESS 


0 . 4740E+05 
0 . 2557E+04 
0 . 2557E+04 
0 . OOOOE+OO 
0 . OOOOE+OO 
0. OOOOE+OO 


0 . 2557E+04 
0 . 62T9E+04 
0 . 1469E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0 . 2557E+04 
0 . 1469E+04 
0 . 6279E+04 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 73T4E+03 
0. OOOOE+OO 
0. OOOOE+OO 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1864E+04 
0. OOOOE+OO 


EFFECTIVE MODULI 


El 1S= 

0 . 4572E+05 

N12S= 

0 . 3300E+00 

E22S= 

0 . 5852E+04 

N23S= 

0 . 2168E+00 

E33S= 

0 . 5852E+04 

G23S= 

0 . 7374E+03 

G13S= 

0 . 1864E+04 

G12S= 

0 . 1864E+04 


0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0. OOOOE+OO 
0 . 1864E+04 
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